random.py revision 3fa19d7ff89be87139e2864fb9186b424d180a58
1e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum"""Random variable generators. 2e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum 3d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters integers 4d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters -------- 5d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters uniform within range 6d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 7d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters sequences 8d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters --------- 9d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters pick random element 10f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger pick random sample 11d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters generate random permutation 12d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 13e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum distributions on the real line: 14e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum ------------------------------ 15d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters uniform 16e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum normal (Gaussian) 17e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum lognormal 18e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum negative exponential 19e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum gamma 20e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum beta 2140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger pareto 2240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Weibull 23e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum 24e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum distributions on the circle (angles 0 to 2pi) 25e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum --------------------------------------------- 26e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum circular uniform 27e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum von Mises 28e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum 2940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond HettingerGeneral notes on the underlying Mersenne Twister core generator: 3040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 3140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger* The period is 2**19937-1. 3240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger* It is one of the most extensively tested generators in existence 3340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger* Without a direct way to compute N steps forward, the 3440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger semantics of jumpahead(n) are weakened to simply jump 3540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger to another distant state and rely on the large period 3640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger to avoid overlapping sequences. 3740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger* The random() method is implemented in C, executes in 3840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a single Python step, and is, therefore, threadsafe. 3940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 40e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum""" 41d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum 422f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettingerfrom warnings import warn as _warn 432f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettingerfrom types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType 44d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersfrom math import log as _log, exp as _exp, pi as _pi, e as _e 45d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersfrom math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin 463fa19d7ff89be87139e2864fb9186b424d180a58Raymond Hettingerfrom math import floor as _floor, ldexp as _ldexp 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", 53356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger "HardwareRandom"] 54ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 55d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersNV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) 56d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersTWOPI = 2.0*_pi 57d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersLOG4 = _log(4.0) 58d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersSG_MAGICCONST = 1.0 + _log(4.5) 592f726e9093381572b21edbfc42659ea89fbdf686Raymond HettingerBPF = 53 # Number of bits in a float 6033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 61356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettingertry: 62356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger from os import urandom as _urandom 63356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger from binascii import hexlify as _hexlify 64356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettingerexcept ImportError: 65356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger _urandom = None 66356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 67356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 68d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters# Translated by Guido van Rossum from C source provided by 6940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# Adrian Baddeley. Adapted by Raymond Hettinger for use with 703fa19d7ff89be87139e2864fb9186b424d180a58Raymond Hettinger# the Mersenne Twister and os.urandom() core generators. 7133d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 72145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettingerimport _random 7340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 74145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettingerclass Random(_random.Random): 75c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Random number generator base class used by bound module functions. 76c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 77c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Used to instantiate instances of Random to get generators that don't 78c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger share state. Especially useful for multi-threaded programs, creating 79c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger a different instance of Random for each thread, and using the jumpahead() 80c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger method to ensure that the generated sequences seen by each thread don't 81c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger overlap. 82c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 83c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Class Random can also be subclassed if you want to use a different basic 84c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger generator of your own devising: in that case, override the following 85c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger methods: random(), seed(), getstate(), setstate() and jumpahead(). 862f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger Optionally, implement a getrandombits() method so that randrange() 872f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger can cover arbitrarily large ranges. 88ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 89c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 9033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 9140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger VERSION = 2 # used by getstate/setstate 9233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 93d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def __init__(self, x=None): 94d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Initialize an instance. 9533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 96d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional argument x controls seeding, as for Random.seed(). 97d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 9833d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 99d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.seed(x) 10040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.gauss_next = None 101ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 1020de88fc4b108751b86443852b6741680d704168fTim Peters def seed(self, a=None): 1030de88fc4b108751b86443852b6741680d704168fTim Peters """Initialize internal state from hashable object. 104d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 105356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger None or no argument seeds from current time or from a hardware 106356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger randomness source if available. 1070de88fc4b108751b86443852b6741680d704168fTim Peters 108bcd725fc456faca13f4598f87c0517f917711cdaTim Peters If a is not None or an int or long, hash(a) is used instead. 109d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 110d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 1113081d59f920229b26293c7a3ee3f1a9da0da53e9Raymond Hettinger if a is None: 112356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger if _urandom is None: 113356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger import time 114356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger a = long(time.time() * 256) # use fractional seconds 115356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger else: 116356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger a = long(_hexlify(_urandom(16)), 16) 117356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 118145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger super(Random, self).seed(a) 11946c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters self.gauss_next = None 12046c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters 121d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def getstate(self): 122d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Return internal state; can be passed to setstate() later.""" 123145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger return self.VERSION, super(Random, self).getstate(), self.gauss_next 124d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 125d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def setstate(self, state): 126d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Restore internal state from object returned by getstate().""" 127d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters version = state[0] 12840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if version == 2: 12940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger version, internalstate, self.gauss_next = state 130145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger super(Random, self).setstate(internalstate) 131d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 132d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError("state with version %s passed to " 133d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters "Random.setstate() of version %s" % 134d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters (version, self.VERSION)) 135d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 136cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- Methods below this point do not need to be overridden when 137cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- subclassing for the purpose of using a different core generator. 138d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 139cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- pickle support ------------------- 140d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 141cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __getstate__(self): # for pickle 142cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters return self.getstate() 143d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 144cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __setstate__(self, state): # for pickle 145cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters self.setstate(state) 146cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 1475f078ff7f0c6bb5086fae077379fc79729c34d2dRaymond Hettinger def __reduce__(self): 1485f078ff7f0c6bb5086fae077379fc79729c34d2dRaymond Hettinger return self.__class__, (), self.getstate() 1495f078ff7f0c6bb5086fae077379fc79729c34d2dRaymond Hettinger 150cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- integer methods ------------------- 151d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 1522f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger def randrange(self, start, stop=None, step=1, int=int, default=None, 1532f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger maxwidth=1L<<BPF): 154d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random item from range(start, stop[, step]). 155d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 156d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters This fixes the problem with randint() which includes the 157d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters endpoint; in Python this is usually not what you want. 1582f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger Do not supply the 'int', 'default', and 'maxwidth' arguments. 159d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 160d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 161d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # This code is a bit messy to make it fast for the 1629146f27b7799dab231083f194a14c6157b57549fTim Peters # common case while still doing adequate error checking. 163d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istart = int(start) 164d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart != start: 165d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer arg 1 for randrange()" 166d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if stop is default: 167d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart > 0: 1682f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if istart >= maxwidth: 1692f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger return self._randbelow(istart) 170d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return int(self.random() * istart) 171d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 1729146f27b7799dab231083f194a14c6157b57549fTim Peters 1739146f27b7799dab231083f194a14c6157b57549fTim Peters # stop argument supplied. 174d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istop = int(stop) 175d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istop != stop: 176d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer stop for randrange()" 1772f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger width = istop - istart 1782f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if step == 1 and width > 0: 17976ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # Note that 1802f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger # int(istart + self.random()*width) 18176ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # instead would be incorrect. For example, consider istart 18276ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # = -2 and istop = 0. Then the guts would be in 18376ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # -2.0 to 0.0 exclusive on both ends (ignoring that random() 18476ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # might return 0.0), and because int() truncates toward 0, the 18576ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # final result would be -1 or 0 (instead of -2 or -1). 1862f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger # istart + int(self.random()*width) 18776ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # would also be incorrect, for a subtler reason: the RHS 18876ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # can return a long, and then randrange() would also return 18976ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # a long, but we're supposed to return an int (for backward 19076ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # compatibility). 1912f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger 1922f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if width >= maxwidth: 19358eb11cf62dd04ccc2c364b62fd51b4265e2e203Tim Peters return int(istart + self._randbelow(width)) 1942f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger return int(istart + int(self.random()*width)) 195d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if step == 1: 1962f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width) 1979146f27b7799dab231083f194a14c6157b57549fTim Peters 1989146f27b7799dab231083f194a14c6157b57549fTim Peters # Non-unit step argument supplied. 199d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istep = int(step) 200d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep != step: 201d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer step for randrange()" 202d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep > 0: 2032f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger n = (width + istep - 1) / istep 204d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif istep < 0: 2052f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger n = (width + istep + 1) / istep 206d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 207d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "zero step for randrange()" 208d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 209d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if n <= 0: 210d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 2112f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger 2122f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if n >= maxwidth: 2132f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger return istart + self._randbelow(n) 214d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return istart + istep*int(self.random() * n) 215d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 216d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def randint(self, a, b): 217cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters """Return random integer in range [a, b], including both end points. 218d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 219d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 220d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return self.randrange(a, b+1) 221d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 2222f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger def _randbelow(self, n, _log=_log, int=int, _maxwidth=1L<<BPF, 2232f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType): 2242f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger """Return a random int in the range [0,n) 2252f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger 2262f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger Handles the case where n has more bits than returned 2272f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger by a single call to the underlying generator. 2282f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger """ 2292f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger 2302f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger try: 2312f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger getrandbits = self.getrandbits 2322f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger except AttributeError: 2332f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger pass 2342f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger else: 2352f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger # Only call self.getrandbits if the original random() builtin method 2362f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger # has not been overridden or if a new getrandbits() was supplied. 2372f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger # This assures that the two methods correspond. 2382f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method: 2392f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger k = int(1.00001 + _log(n-1, 2.0)) # 2**k > n-1 > 2**(k-2) 2402f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger r = getrandbits(k) 2412f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger while r >= n: 2422f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger r = getrandbits(k) 2432f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger return r 2442f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if n >= _maxwidth: 2452f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger _warn("Underlying random() generator does not supply \n" 2462f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger "enough bits to choose from a population range this large") 2472f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger return int(self.random() * n) 2482f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger 249cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- sequence methods ------------------- 250cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 251d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def choice(self, seq): 252d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random element from a non-empty sequence.""" 2535dae505bbd59641a948c81bea981e7c44d4c2343Raymond Hettinger return seq[int(self.random() * len(seq))] # raises IndexError if seq is empty 254d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 255d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def shuffle(self, x, random=None, int=int): 256d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """x, random=random.random -> shuffle list x in place; return None. 257d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 258d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional arg random is a 0-argument function returning a random 259d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters float in [0.0, 1.0); by default, the standard random.random. 260d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 261d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Note that for even rather small len(x), the total number of 262d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters permutations of x is larger than the period of most random number 263d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters generators; this implies that "most" permutations of a long 264d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters sequence can never be generated. 265d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 266d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 267d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if random is None: 268d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 26985c20a41dfcec04d161ad7da7260e7b94c62d228Raymond Hettinger for i in reversed(xrange(1, len(x))): 270cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # pick an element in x[:i+1] with which to exchange x[i] 271d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters j = int(random() * (i+1)) 272d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x[i], x[j] = x[j], x[i] 273d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 274fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger def sample(self, population, k): 275f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger """Chooses k unique random elements from a population sequence. 276f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 277c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger Returns a new list containing elements from the population while 278c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger leaving the original population unchanged. The resulting list is 279c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger in selection order so that all sub-slices will also be valid random 280c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger samples. This allows raffle winners (the sample) to be partitioned 281c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger into grand prize and second place winners (the subslices). 282f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 283c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger Members of the population need not be hashable or unique. If the 284c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger population contains repeats, then each occurrence is a possible 285c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger selection in the sample. 286f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 287c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger To choose a sample in a range of integers, use xrange as an argument. 288c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger This is especially fast and space efficient for sampling from a 289c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger large population: sample(xrange(10000000), 60) 290f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger """ 291f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 292c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger # Sampling without replacement entails tracking either potential 2938b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # selections (the pool) in a list or previous selections in a 2948b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # dictionary. 295c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger 2962b55d35850e3e8e0b28aba7878d3f9122a7907acJeremy Hylton # When the number of selections is small compared to the 2972b55d35850e3e8e0b28aba7878d3f9122a7907acJeremy Hylton # population, then tracking selections is efficient, requiring 2982b55d35850e3e8e0b28aba7878d3f9122a7907acJeremy Hylton # only a small dictionary and an occasional reselection. For 2992b55d35850e3e8e0b28aba7878d3f9122a7907acJeremy Hylton # a larger number of selections, the pool tracking method is 3002b55d35850e3e8e0b28aba7878d3f9122a7907acJeremy Hylton # preferred since the list takes less space than the 3012b55d35850e3e8e0b28aba7878d3f9122a7907acJeremy Hylton # dictionary and it doesn't suffer from frequent reselections. 302c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger 303f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger n = len(population) 304f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger if not 0 <= k <= n: 305f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger raise ValueError, "sample larger than population" 3068b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger random = self.random 307fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger _int = int 308c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger result = [None] * k 309f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger if n < 6 * k: # if n len list takes less space than a k len dict 310311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger pool = list(population) 311311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger for i in xrange(k): # invariant: non-selected at [0,n-i) 312fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger j = _int(random() * (n-i)) 313311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger result[i] = pool[j] 3148b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger pool[j] = pool[n-i-1] # move non-selected item into vacancy 315c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger else: 31666d09f1b3029d9cf975ccf26c437c9fb2605db91Raymond Hettinger try: 31766d09f1b3029d9cf975ccf26c437c9fb2605db91Raymond Hettinger n > 0 and (population[0], population[n//2], population[n-1]) 31866d09f1b3029d9cf975ccf26c437c9fb2605db91Raymond Hettinger except (TypeError, KeyError): # handle sets and dictionaries 31966d09f1b3029d9cf975ccf26c437c9fb2605db91Raymond Hettinger population = tuple(population) 320311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger selected = {} 321c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger for i in xrange(k): 322fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger j = _int(random() * n) 323311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while j in selected: 324fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger j = _int(random() * n) 325c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger result[i] = selected[j] = population[j] 326311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger return result 327f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 328cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- real-valued distributions ------------------- 329cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 330cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- uniform distribution ------------------- 331d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 332d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def uniform(self, a, b): 333d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Get a random number in the range [a, b).""" 334d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return a + (b-a) * self.random() 335ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 336cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- normal distribution -------------------- 337ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 338d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def normalvariate(self, mu, sigma): 339c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Normal distribution. 340c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 341c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean, and sigma is the standard deviation. 342ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 343c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 344d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mu = mean, sigma = standard deviation 345d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 346d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses Kinderman and Monahan method. Reference: Kinderman, 347d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # A.J. and Monahan, J.F., "Computer generation of random 348d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables using the ratio of uniform deviates", ACM Trans 349d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Math Software, 3, (1977), pp257-260. 350d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 351d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 352311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while True: 353d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 35473ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u2 = 1.0 - random() 355d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = NV_MAGICCONST*(u1-0.5)/u2 356d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters zz = z*z/4.0 357d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if zz <= -_log(u2): 358d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 359d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 360ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 361cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- lognormal distribution -------------------- 362ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 363d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def lognormvariate(self, mu, sigma): 364c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Log normal distribution. 365c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 366c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger If you take the natural logarithm of this distribution, you'll get a 367c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger normal distribution with mean mu and standard deviation sigma. 368c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu can have any value, and sigma must be greater than zero. 369ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 370c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 371d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return _exp(self.normalvariate(mu, sigma)) 372ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 373cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- exponential distribution -------------------- 374ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 375d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def expovariate(self, lambd): 376c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Exponential distribution. 377c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 378c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger lambd is 1.0 divided by the desired mean. (The parameter would be 379c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger called "lambda", but that is a reserved word in Python.) Returned 380c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger values range from 0 to positive infinity. 381ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 382c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 383d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lambd: rate lambd = 1/mean 384d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # ('lambda' is a Python reserved word) 385ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 386d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 3870c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 388d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 389d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 390d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return -_log(u)/lambd 391ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 392cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- von Mises distribution -------------------- 393ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 394d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def vonmisesvariate(self, mu, kappa): 395c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Circular data distribution. 396ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 397c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean angle, expressed in radians between 0 and 2*pi, and 398c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger kappa is the concentration parameter, which must be greater than or 399c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger equal to zero. If kappa is equal to zero, this distribution reduces 400c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger to a uniform random angle over the range 0 to 2*pi. 401ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 402c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 403d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mu: mean angle (in radians between 0 and 2*pi) 404d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # kappa: concentration parameter kappa (>= 0) 405d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # if kappa = 0 generate uniform random angle 4065810297052003f28788f6790ac799fe8e5373494Guido van Rossum 407d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Based upon an algorithm published in: Fisher, N.I., 408d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # "Statistical Analysis of Circular Data", Cambridge 409d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # University Press, 1993. 4105810297052003f28788f6790ac799fe8e5373494Guido van Rossum 411d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Thanks to Magnus Kessler for a correction to the 412d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # implementation of step 4. 4135810297052003f28788f6790ac799fe8e5373494Guido van Rossum 414d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 415d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if kappa <= 1e-6: 416d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return TWOPI * random() 417ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 418d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa) 419d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (a - _sqrt(2.0 * a))/(2.0 * kappa) 420d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = (1.0 + b * b)/(2.0 * b) 421ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 422311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while True: 423d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 424ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 425d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(_pi * u1) 426d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters f = (1.0 + r * z)/(r + z) 427d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters c = kappa * (r - f) 428ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 429d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u2 = random() 430ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 431d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)): 432d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 433ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 434d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u3 = random() 435d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if u3 > 0.5: 436d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters theta = (mu % TWOPI) + _acos(f) 437d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 438d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters theta = (mu % TWOPI) - _acos(f) 439ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 440d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return theta 441ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 442cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- gamma distribution -------------------- 443ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 444d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gammavariate(self, alpha, beta): 445c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Gamma distribution. Not the gamma function! 446c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 447c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Conditions on the parameters are alpha > 0 and beta > 0. 448c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 449c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 4508ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 451b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 4528ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 453570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum # Warning: a few older sources define the gamma distribution in terms 454570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum # of alpha > -1.0 455570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum if alpha <= 0.0 or beta <= 0.0: 456570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum raise ValueError, 'gammavariate: alpha and beta must be > 0.0' 4578ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 458d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 459d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if alpha > 1.0: 460d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 461d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses R.C.H. Cheng, "The generation of Gamma 462d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables with non-integral shape parameters", 463d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Applied Statistics, (1977), 26, No. 1, p71-74 464d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 465ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger ainv = _sqrt(2.0 * alpha - 1.0) 466ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger bbb = alpha - LOG4 467ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger ccc = alpha + ainv 4688ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 469311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while True: 470d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 47173ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger if not 1e-7 < u1 < .9999999: 47273ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger continue 47373ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u2 = 1.0 - random() 474d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters v = _log(u1/(1.0-u1))/ainv 475d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = alpha*_exp(v) 476d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = u1*u1*u2 477d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = bbb+ccc*v-x 478d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 479b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return x * beta 480d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 481d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif alpha == 1.0: 482d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # expovariate(1) 4830c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 484d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 485d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 486b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return -_log(u) * beta 487d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 488d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: # alpha is between 0 and 1 (exclusive) 489d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 490d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 491d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 492311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while True: 493d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 494d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (_e + alpha)/_e 495d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters p = b*u 496d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if p <= 1.0: 497d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = pow(p, 1.0/alpha) 498d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 499d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # p > 1 500d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = -_log((b-p)/alpha) 501d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 502d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not (((p <= 1.0) and (u1 > _exp(-x))) or 503d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters ((p > 1) and (u1 > pow(x, alpha - 1.0)))): 504d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 505b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return x * beta 506b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger 507cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Gauss (faster alternative) -------------------- 50895bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 509d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gauss(self, mu, sigma): 510c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Gaussian distribution. 511c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 512c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean, and sigma is the standard deviation. This is 513c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger slightly faster than the normalvariate() function. 514c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 515c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Not thread-safe without a lock around calls. 516ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 517c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 518d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 519d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # When x and y are two variables from [0, 1), uniformly 520d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # distributed, then 521d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 522d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # cos(2*pi*x)*sqrt(-2*log(1-y)) 523d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # sin(2*pi*x)*sqrt(-2*log(1-y)) 524d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 525d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # are two *independent* variables with normal distribution 526d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (mu = 0, sigma = 1). 527d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (Lambert Meertens) 528d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (corrected version; bug discovered by Mike Miller, fixed by LM) 529d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 530d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Multithreading note: When two threads call this function 531d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # simultaneously, it is possible that they will receive the 532d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # same return value. The window is very small though. To 533d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # avoid this, you have to use a lock around all calls. (I 534d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # didn't want to slow this down in the serial case by using a 535d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lock here.) 536d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 537d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 538d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = self.gauss_next 539d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = None 540d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if z is None: 541d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x2pi = random() * TWOPI 542d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters g2rad = _sqrt(-2.0 * _log(1.0 - random())) 543d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(x2pi) * g2rad 544d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = _sin(x2pi) * g2rad 545d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 546d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 54795bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 548cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- beta -------------------- 54985e2e4742d0a1accecd02058a7907df36308297eTim Peters## See 55085e2e4742d0a1accecd02058a7907df36308297eTim Peters## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470 55185e2e4742d0a1accecd02058a7907df36308297eTim Peters## for Ivan Frohne's insightful analysis of why the original implementation: 55285e2e4742d0a1accecd02058a7907df36308297eTim Peters## 55385e2e4742d0a1accecd02058a7907df36308297eTim Peters## def betavariate(self, alpha, beta): 55485e2e4742d0a1accecd02058a7907df36308297eTim Peters## # Discrete Event Simulation in C, pp 87-88. 55585e2e4742d0a1accecd02058a7907df36308297eTim Peters## 55685e2e4742d0a1accecd02058a7907df36308297eTim Peters## y = self.expovariate(alpha) 55785e2e4742d0a1accecd02058a7907df36308297eTim Peters## z = self.expovariate(1.0/beta) 55885e2e4742d0a1accecd02058a7907df36308297eTim Peters## return z/(y+z) 55985e2e4742d0a1accecd02058a7907df36308297eTim Peters## 56085e2e4742d0a1accecd02058a7907df36308297eTim Peters## was dead wrong, and how it probably got that way. 56195bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 562d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def betavariate(self, alpha, beta): 563c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Beta distribution. 564c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 565c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Conditions on the parameters are alpha > -1 and beta} > -1. 566c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Returned values range between 0 and 1. 567ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 568c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 569ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 57085e2e4742d0a1accecd02058a7907df36308297eTim Peters # This version due to Janne Sinkkonen, and matches all the std 57185e2e4742d0a1accecd02058a7907df36308297eTim Peters # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 57285e2e4742d0a1accecd02058a7907df36308297eTim Peters y = self.gammavariate(alpha, 1.) 57385e2e4742d0a1accecd02058a7907df36308297eTim Peters if y == 0: 57485e2e4742d0a1accecd02058a7907df36308297eTim Peters return 0.0 57585e2e4742d0a1accecd02058a7907df36308297eTim Peters else: 57685e2e4742d0a1accecd02058a7907df36308297eTim Peters return y / (y + self.gammavariate(beta, 1.)) 57795bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 578cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Pareto -------------------- 579cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 580d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def paretovariate(self, alpha): 581c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Pareto distribution. alpha is the shape parameter.""" 582d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 495 583cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 58473ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u = 1.0 - self.random() 585d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return 1.0 / pow(u, 1.0/alpha) 586cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 587cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Weibull -------------------- 588cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 589d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def weibullvariate(self, alpha, beta): 590c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Weibull distribution. 591c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 592c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger alpha is the scale parameter and beta is the shape parameter. 593ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 594c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 595d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 499; bug fix courtesy Bill Arms 596cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 59773ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u = 1.0 - self.random() 598d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return alpha * pow(-_log(u), 1.0/beta) 5996c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum 60040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger## -------------------- Wichmann-Hill ------------------- 60140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 60240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettingerclass WichmannHill(Random): 60340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 60440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger VERSION = 1 # used by getstate/setstate 60540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 60640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def seed(self, a=None): 60740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Initialize internal state from hashable object. 60840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 609356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger None or no argument seeds from current time or from a hardware 610356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger randomness source if available. 61140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 61240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger If a is not None or an int or long, hash(a) is used instead. 61340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 61440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger If a is an int or long, a is used directly. Distinct values between 61540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 0 and 27814431486575L inclusive are guaranteed to yield distinct 61640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger internal states (this guarantee is specific to the default 61740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Wichmann-Hill generator). 61840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 61940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 62040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if a is None: 621356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger if _urandom is None: 622356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger import time 623356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger a = long(time.time() * 256) # use fractional seconds 624356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger else: 625356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger a = long(_hexlify(_urandom(16)), 16) 62640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 62740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not isinstance(a, (int, long)): 62840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a = hash(a) 62940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 63040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, x = divmod(a, 30268) 63140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, y = divmod(a, 30306) 63240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, z = divmod(a, 30322) 63340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = int(x)+1, int(y)+1, int(z)+1 63440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 63540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.gauss_next = None 63640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 63740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def random(self): 63840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Get the next random number in the range [0.0, 1.0).""" 63940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 64040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Wichman-Hill random number generator. 64140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # 64240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Wichmann, B. A. & Hill, I. D. (1982) 64340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Algorithm AS 183: 64440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # An efficient and portable pseudo-random number generator 64540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Applied Statistics 31 (1982) 188-190 64640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # 64740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # see also: 64840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Correction to Algorithm AS 183 64940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Applied Statistics 33 (1984) 123 65040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # 65140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # McLeod, A. I. (1985) 65240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # A remark on Algorithm AS 183 65340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Applied Statistics 34 (1985),198-200 65440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 65540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # This part is thread-unsafe: 65640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # BEGIN CRITICAL SECTION 65740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x, y, z = self._seed 65840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x = (171 * x) % 30269 65940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger y = (172 * y) % 30307 66040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger z = (170 * z) % 30323 66140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = x, y, z 66240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # END CRITICAL SECTION 66340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 66440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Note: on a platform using IEEE-754 double arithmetic, this can 66540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # never return 0.0 (asserted by Tim; proof too long for a comment). 66640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0 66740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 66840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def getstate(self): 66940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Return internal state; can be passed to setstate() later.""" 67040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger return self.VERSION, self._seed, self.gauss_next 67140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 67240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def setstate(self, state): 67340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Restore internal state from object returned by getstate().""" 67440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger version = state[0] 67540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if version == 1: 67640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger version, self._seed, self.gauss_next = state 67740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger else: 67840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise ValueError("state with version %s passed to " 67940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger "Random.setstate() of version %s" % 68040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger (version, self.VERSION)) 68140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 68240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def jumpahead(self, n): 68340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Act as if n calls to random() were made, but quickly. 68440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 68540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger n is an int, greater than or equal to 0. 68640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 68740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Example use: If you have 2 threads and know that each will 68840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger consume no more than a million random numbers, create two Random 68940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger objects r1 and r2, then do 69040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger r2.setstate(r1.getstate()) 69140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger r2.jumpahead(1000000) 69240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Then r1 and r2 will use guaranteed-disjoint segments of the full 69340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger period. 69440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 69540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 69640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not n >= 0: 69740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise ValueError("n must be >= 0") 69840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x, y, z = self._seed 69940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x = int(x * pow(171, n, 30269)) % 30269 70040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger y = int(y * pow(172, n, 30307)) % 30307 70140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger z = int(z * pow(170, n, 30323)) % 30323 70240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = x, y, z 70340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 70440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def __whseed(self, x=0, y=0, z=0): 70540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Set the Wichmann-Hill seed from (x, y, z). 70640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 70740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger These must be integers in the range [0, 256). 70840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 70940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 71040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not type(x) == type(y) == type(z) == int: 71140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise TypeError('seeds must be integers') 71240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): 71340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise ValueError('seeds must be in range(0, 256)') 71440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if 0 == x == y == z: 71540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Initialize from current time 71640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger import time 71740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t = long(time.time() * 256) 71840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t = int((t&0xffffff) ^ (t>>24)) 71940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t, x = divmod(t, 256) 72040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t, y = divmod(t, 256) 72140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t, z = divmod(t, 256) 72240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Zero is a poor seed, so substitute 1 72340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = (x or 1, y or 1, z or 1) 72440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 72540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.gauss_next = None 72640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 72740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def whseed(self, a=None): 72840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Seed from hashable object's hash code. 72940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 73040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger None or no argument seeds from current time. It is not guaranteed 73140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger that objects with distinct hash codes lead to distinct internal 73240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger states. 73340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 73440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger This is obsolete, provided for compatibility with the seed routine 73540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger used prior to Python 2.1. Use the .seed() method instead. 73640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 73740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 73840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if a is None: 73940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.__whseed() 74040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger return 74140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a = hash(a) 74240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, x = divmod(a, 256) 74340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, y = divmod(a, 256) 74440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, z = divmod(a, 256) 74540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x = (x + a) % 256 or 1 74640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger y = (y + a) % 256 or 1 74740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger z = (z + a) % 256 or 1 74840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.__whseed(x, y, z) 74940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 750356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger## -------------------- Hardware Random Source ------------------- 751356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 752356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettingerclass HardwareRandom(Random): 753356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger """Alternate random number generator using hardware sources. 754356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 755356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger Not available on all systems (see os.urandom() for details). 756356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger """ 757356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 758356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger def random(self): 759356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger """Get the next random number in the range [0.0, 1.0).""" 760356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger if _urandom is None: 761356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger raise NotImplementedError('Cannot find hardware entropy source') 7623fa19d7ff89be87139e2864fb9186b424d180a58Raymond Hettinger return _ldexp(long(_hexlify(_urandom(7)), 16) >> 3, -BPF) 763356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 764356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger def getrandbits(self, k): 765356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger """getrandbits(k) -> x. Generates a long int with k random bits.""" 766356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger if _urandom is None: 767356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger raise NotImplementedError('Cannot find hardware entropy source') 768356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger if k <= 0: 769356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger raise ValueError('number of bits must be greater than zero') 770356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger if k != int(k): 771356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger raise TypeError('number of bits should be an integer') 772356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger bytes = (k + 7) // 8 # bits / 8 and rounded up 773356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger x = long(_hexlify(_urandom(bytes)), 16) 774356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger return x >> (bytes * 8 - k) # trim excess bits 775356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 776356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger def _stub(self, *args, **kwds): 777356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger "Stub method. Not used for a hardware random number generator." 778356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger return None 779356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger seed = jumpahead = _stub 780356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 781356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger def _notimplemented(self, *args, **kwds): 782356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger "Method should not be called for a hardware random number generator." 783356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger raise NotImplementedError('Hardware entropy source does not have state.') 784356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger getstate = setstate = _notimplemented 785356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 786cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- test program -------------------- 787ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 78862297132215490e9cb406e1a21f03aff40d421cbRaymond Hettingerdef _test_generator(n, func, args): 7890c9886d589ddebf32de0ca3f027a173222ed383aTim Peters import time 79062297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger print n, 'times', func.__name__ 791b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger total = 0.0 7920c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = 0.0 7930c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = 1e10 7940c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = -1e10 7950c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t0 = time.time() 7960c9886d589ddebf32de0ca3f027a173222ed383aTim Peters for i in range(n): 79762297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger x = func(*args) 798b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger total += x 7990c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = sqsum + x*x 8000c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = min(x, smallest) 8010c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = max(x, largest) 8020c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t1 = time.time() 8030c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print round(t1-t0, 3), 'sec,', 804b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger avg = total/n 805d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters stddev = _sqrt(sqsum/n - avg*avg) 8060c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print 'avg %g, stddev %g, min %g, max %g' % \ 8070c9886d589ddebf32de0ca3f027a173222ed383aTim Peters (avg, stddev, smallest, largest) 808ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 809f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 810f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettingerdef _test(N=2000): 81162297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, random, ()) 81262297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, normalvariate, (0.0, 1.0)) 81362297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, lognormvariate, (0.0, 1.0)) 81462297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, vonmisesvariate, (0.0, 1.0)) 81562297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (0.01, 1.0)) 81662297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (0.1, 1.0)) 81762297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (0.1, 2.0)) 81862297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (0.5, 1.0)) 81962297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (0.9, 1.0)) 82062297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (1.0, 1.0)) 82162297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (2.0, 1.0)) 82262297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (20.0, 1.0)) 82362297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (200.0, 1.0)) 82462297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gauss, (0.0, 1.0)) 82562297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, betavariate, (3.0, 3.0)) 826cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 827715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# Create one instance, seeded from current time, and export its methods 82840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# as module-level functions. The functions share state across all uses 82940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger#(both in the user's code and in the Python libraries), but that's fine 83040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# for most programs and is easier for the casual user than making them 83140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# instantiate their own Random() instance. 83240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 833d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_inst = Random() 834d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersseed = _inst.seed 835d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandom = _inst.random 836d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersuniform = _inst.uniform 837d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandint = _inst.randint 838d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterschoice = _inst.choice 839d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandrange = _inst.randrange 840f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettingersample = _inst.sample 841d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersshuffle = _inst.shuffle 842d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersnormalvariate = _inst.normalvariate 843d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterslognormvariate = _inst.lognormvariate 844d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersexpovariate = _inst.expovariate 845d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersvonmisesvariate = _inst.vonmisesvariate 846d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgammavariate = _inst.gammavariate 847d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgauss = _inst.gauss 848d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersbetavariate = _inst.betavariate 849d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersparetovariate = _inst.paretovariate 850d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersweibullvariate = _inst.weibullvariate 851d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgetstate = _inst.getstate 852d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterssetstate = _inst.setstate 853d52269bfd029c4a517ea74c17edd5c3a250c366cTim Petersjumpahead = _inst.jumpahead 8542f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettingergetrandbits = _inst.getrandbits 855d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 856ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumif __name__ == '__main__': 857d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test() 858