14adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao"""Random variable generators. 24adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 34adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao integers 44adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao -------- 54adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao uniform within range 64adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 74adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao sequences 84adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao --------- 94adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao pick random element 104adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao pick random sample 114adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao generate random permutation 124adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 134adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao distributions on the real line: 144adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao ------------------------------ 154adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao uniform 164adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao triangular 174adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao normal (Gaussian) 184adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao lognormal 194adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao negative exponential 204adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao gamma 214adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao beta 224adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao pareto 234adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Weibull 244adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 254adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao distributions on the circle (angles 0 to 2pi) 264adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao --------------------------------------------- 274adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao circular uniform 284adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao von Mises 294adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 304adfde8bc82dd39f59e0445588c3e599ada477dJosh GaoGeneral notes on the underlying Mersenne Twister core generator: 314adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 324adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao* The period is 2**19937-1. 334adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao* It is one of the most extensively tested generators in existence. 344adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao* Without a direct way to compute N steps forward, the semantics of 354adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao jumpahead(n) are weakened to simply jump to another distant state and rely 364adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao on the large period to avoid overlapping sequences. 374adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao* The random() method is implemented in C, executes in a single Python step, 384adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao and is, therefore, threadsafe. 394adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 404adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao""" 414adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 424adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaofrom __future__ import division 434adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaofrom warnings import warn as _warn 444adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaofrom types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType 454adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaofrom math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil 464adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaofrom math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin 474adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaofrom os import urandom as _urandom 484adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaofrom binascii import hexlify as _hexlify 494adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaoimport hashlib as _hashlib 504adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 514adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao__all__ = ["Random","seed","random","uniform","randint","choice","sample", 524adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao "randrange","shuffle","normalvariate","lognormvariate", 534adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao "expovariate","vonmisesvariate","gammavariate","triangular", 544adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao "gauss","betavariate","paretovariate","weibullvariate", 554adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao "getstate","setstate","jumpahead", "WichmannHill", "getrandbits", 564adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao "SystemRandom"] 574adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 584adfde8bc82dd39f59e0445588c3e599ada477dJosh GaoNV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) 594adfde8bc82dd39f59e0445588c3e599ada477dJosh GaoTWOPI = 2.0*_pi 604adfde8bc82dd39f59e0445588c3e599ada477dJosh GaoLOG4 = _log(4.0) 614adfde8bc82dd39f59e0445588c3e599ada477dJosh GaoSG_MAGICCONST = 1.0 + _log(4.5) 624adfde8bc82dd39f59e0445588c3e599ada477dJosh GaoBPF = 53 # Number of bits in a float 634adfde8bc82dd39f59e0445588c3e599ada477dJosh GaoRECIP_BPF = 2**-BPF 644adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 654adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 664adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao# Translated by Guido van Rossum from C source provided by 674adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao# Adrian Baddeley. Adapted by Raymond Hettinger for use with 684adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao# the Mersenne Twister and os.urandom() core generators. 694adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 704adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaoimport _random 714adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 724adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaoclass Random(_random.Random): 734adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Random number generator base class used by bound module functions. 744adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 754adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Used to instantiate instances of Random to get generators that don't 764adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao share state. Especially useful for multi-threaded programs, creating 774adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao a different instance of Random for each thread, and using the jumpahead() 784adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao method to ensure that the generated sequences seen by each thread don't 794adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao overlap. 804adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 814adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Class Random can also be subclassed if you want to use a different basic 824adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao generator of your own devising: in that case, override the following 834adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao methods: random(), seed(), getstate(), setstate() and jumpahead(). 844adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Optionally, implement a getrandbits() method so that randrange() can cover 854adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao arbitrarily large ranges. 864adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 874adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 884adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 894adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao VERSION = 3 # used by getstate/setstate 904adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 914adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def __init__(self, x=None): 924adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Initialize an instance. 934adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 944adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Optional argument x controls seeding, as for Random.seed(). 954adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 964adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 974adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao self.seed(x) 984adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao self.gauss_next = None 994adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1004adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def seed(self, a=None): 1014adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Initialize internal state from hashable object. 1024adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1034adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao None or no argument seeds from current time or from an operating 1044adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao system specific randomness source if available. 1054adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1064adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao If a is not None or an int or long, hash(a) is used instead. 1074adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 1084adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1094adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if a is None: 1104adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao try: 1114adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao a = long(_hexlify(_urandom(16)), 16) 1124adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao except NotImplementedError: 1134adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao import time 1144adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao a = long(time.time() * 256) # use fractional seconds 1154adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1164adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao super(Random, self).seed(a) 1174adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao self.gauss_next = None 1184adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1194adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def getstate(self): 1204adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Return internal state; can be passed to setstate() later.""" 1214adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return self.VERSION, super(Random, self).getstate(), self.gauss_next 1224adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1234adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def setstate(self, state): 1244adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Restore internal state from object returned by getstate().""" 1254adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao version = state[0] 1264adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if version == 3: 1274adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao version, internalstate, self.gauss_next = state 1284adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao super(Random, self).setstate(internalstate) 1294adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao elif version == 2: 1304adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao version, internalstate, self.gauss_next = state 1314adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # In version 2, the state was saved as signed ints, which causes 1324adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # inconsistencies between 32/64-bit systems. The state is 1334adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # really unsigned 32-bit ints, so we convert negative ints from 1344adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # version 2 to positive longs for version 3. 1354adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao try: 1364adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao internalstate = tuple( long(x) % (2**32) for x in internalstate ) 1374adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao except ValueError, e: 1384adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise TypeError, e 1394adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao super(Random, self).setstate(internalstate) 1404adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao else: 1414adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise ValueError("state with version %s passed to " 1424adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao "Random.setstate() of version %s" % 1434adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao (version, self.VERSION)) 1444adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1454adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def jumpahead(self, n): 1464adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Change the internal state to one that is likely far away 1474adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao from the current state. This method will not be in Py3.x, 1484adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao so it is better to simply reseed. 1494adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 1504adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # The super.jumpahead() method uses shuffling to change state, 1514adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # so it needs a large and "interesting" n to work with. Here, 1524adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # we use hashing to create a large n for the shuffle. 1534adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao s = repr(n) + repr(self.getstate()) 1544adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao n = int(_hashlib.new('sha512', s).hexdigest(), 16) 1554adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao super(Random, self).jumpahead(n) 1564adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1574adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## ---- Methods below this point do not need to be overridden when 1584adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## ---- subclassing for the purpose of using a different core generator. 1594adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1604adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- pickle support ------------------- 1614adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1624adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def __getstate__(self): # for pickle 1634adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return self.getstate() 1644adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1654adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def __setstate__(self, state): # for pickle 1664adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao self.setstate(state) 1674adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1684adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def __reduce__(self): 1694adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return self.__class__, (), self.getstate() 1704adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1714adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- integer methods ------------------- 1724adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1734adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def randrange(self, start, stop=None, step=1, int=int, default=None, 1744adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao maxwidth=1L<<BPF): 1754adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Choose a random item from range(start, stop[, step]). 1764adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1774adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao This fixes the problem with randint() which includes the 1784adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao endpoint; in Python this is usually not what you want. 1794adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Do not supply the 'int', 'default', and 'maxwidth' arguments. 1804adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 1814adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1824adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # This code is a bit messy to make it fast for the 1834adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # common case while still doing adequate error checking. 1844adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao istart = int(start) 1854adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if istart != start: 1864adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise ValueError, "non-integer arg 1 for randrange()" 1874adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if stop is default: 1884adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if istart > 0: 1894adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if istart >= maxwidth: 1904adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return self._randbelow(istart) 1914adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return int(self.random() * istart) 1924adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise ValueError, "empty range for randrange()" 1934adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 1944adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # stop argument supplied. 1954adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao istop = int(stop) 1964adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if istop != stop: 1974adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise ValueError, "non-integer stop for randrange()" 1984adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao width = istop - istart 1994adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if step == 1 and width > 0: 2004adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Note that 2014adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # int(istart + self.random()*width) 2024adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # instead would be incorrect. For example, consider istart 2034adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # = -2 and istop = 0. Then the guts would be in 2044adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # -2.0 to 0.0 exclusive on both ends (ignoring that random() 2054adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # might return 0.0), and because int() truncates toward 0, the 2064adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # final result would be -1 or 0 (instead of -2 or -1). 2074adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # istart + int(self.random()*width) 2084adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # would also be incorrect, for a subtler reason: the RHS 2094adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # can return a long, and then randrange() would also return 2104adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # a long, but we're supposed to return an int (for backward 2114adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # compatibility). 2124adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2134adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if width >= maxwidth: 2144adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return int(istart + self._randbelow(width)) 2154adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return int(istart + int(self.random()*width)) 2164adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if step == 1: 2174adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width) 2184adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2194adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Non-unit step argument supplied. 2204adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao istep = int(step) 2214adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if istep != step: 2224adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise ValueError, "non-integer step for randrange()" 2234adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if istep > 0: 2244adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao n = (width + istep - 1) // istep 2254adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao elif istep < 0: 2264adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao n = (width + istep + 1) // istep 2274adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao else: 2284adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise ValueError, "zero step for randrange()" 2294adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2304adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if n <= 0: 2314adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise ValueError, "empty range for randrange()" 2324adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2334adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if n >= maxwidth: 2344adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return istart + istep*self._randbelow(n) 2354adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return istart + istep*int(self.random() * n) 2364adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2374adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def randint(self, a, b): 2384adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Return random integer in range [a, b], including both end points. 2394adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 2404adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2414adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return self.randrange(a, b+1) 2424adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2434adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def _randbelow(self, n, _log=_log, int=int, _maxwidth=1L<<BPF, 2444adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType): 2454adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Return a random int in the range [0,n) 2464adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2474adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Handles the case where n has more bits than returned 2484adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao by a single call to the underlying generator. 2494adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 2504adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2514adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao try: 2524adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao getrandbits = self.getrandbits 2534adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao except AttributeError: 2544adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao pass 2554adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao else: 2564adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Only call self.getrandbits if the original random() builtin method 2574adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # has not been overridden or if a new getrandbits() was supplied. 2584adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # This assures that the two methods correspond. 2594adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method: 2604adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao k = int(1.00001 + _log(n-1, 2.0)) # 2**k > n-1 > 2**(k-2) 2614adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao r = getrandbits(k) 2624adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao while r >= n: 2634adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao r = getrandbits(k) 2644adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return r 2654adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if n >= _maxwidth: 2664adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _warn("Underlying random() generator does not supply \n" 2674adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao "enough bits to choose from a population range this large") 2684adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return int(self.random() * n) 2694adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2704adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- sequence methods ------------------- 2714adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2724adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def choice(self, seq): 2734adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Choose a random element from a non-empty sequence.""" 2744adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return seq[int(self.random() * len(seq))] # raises IndexError if seq is empty 2754adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2764adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def shuffle(self, x, random=None, int=int): 2774adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """x, random=random.random -> shuffle list x in place; return None. 2784adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2794adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Optional arg random is a 0-argument function returning a random 2804adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao float in [0.0, 1.0); by default, the standard random.random. 2814adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 2824adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2834adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if random is None: 2844adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao random = self.random 2854adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao for i in reversed(xrange(1, len(x))): 2864adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # pick an element in x[:i+1] with which to exchange x[i] 2874adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao j = int(random() * (i+1)) 2884adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao x[i], x[j] = x[j], x[i] 2894adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2904adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def sample(self, population, k): 2914adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Chooses k unique random elements from a population sequence. 2924adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2934adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Returns a new list containing elements from the population while 2944adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao leaving the original population unchanged. The resulting list is 2954adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao in selection order so that all sub-slices will also be valid random 2964adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao samples. This allows raffle winners (the sample) to be partitioned 2974adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao into grand prize and second place winners (the subslices). 2984adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 2994adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Members of the population need not be hashable or unique. If the 3004adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao population contains repeats, then each occurrence is a possible 3014adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao selection in the sample. 3024adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3034adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao To choose a sample in a range of integers, use xrange as an argument. 3044adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao This is especially fast and space efficient for sampling from a 3054adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao large population: sample(xrange(10000000), 60) 3064adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 3074adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3084adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Sampling without replacement entails tracking either potential 3094adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # selections (the pool) in a list or previous selections in a set. 3104adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3114adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # When the number of selections is small compared to the 3124adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # population, then tracking selections is efficient, requiring 3134adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # only a small set and an occasional reselection. For 3144adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # a larger number of selections, the pool tracking method is 3154adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # preferred since the list takes less space than the 3164adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # set and it doesn't suffer from frequent reselections. 3174adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3184adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao n = len(population) 3194adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if not 0 <= k <= n: 3204adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise ValueError("sample larger than population") 3214adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao random = self.random 3224adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _int = int 3234adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao result = [None] * k 3244adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao setsize = 21 # size of a small set minus size of an empty list 3254adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if k > 5: 3264adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets 3274adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if n <= setsize or hasattr(population, "keys"): 3284adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # An n-length list is smaller than a k-length set, or this is a 3294adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # mapping type so the other algorithm wouldn't work. 3304adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao pool = list(population) 3314adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao for i in xrange(k): # invariant: non-selected at [0,n-i) 3324adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao j = _int(random() * (n-i)) 3334adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao result[i] = pool[j] 3344adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao pool[j] = pool[n-i-1] # move non-selected item into vacancy 3354adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao else: 3364adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao try: 3374adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao selected = set() 3384adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao selected_add = selected.add 3394adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao for i in xrange(k): 3404adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao j = _int(random() * n) 3414adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao while j in selected: 3424adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao j = _int(random() * n) 3434adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao selected_add(j) 3444adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao result[i] = population[j] 3454adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao except (TypeError, KeyError): # handle (at least) sets 3464adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if isinstance(population, list): 3474adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise 3484adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return self.sample(tuple(population), k) 3494adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return result 3504adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3514adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- real-valued distributions ------------------- 3524adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3534adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- uniform distribution ------------------- 3544adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3554adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def uniform(self, a, b): 3564adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao "Get a random number in the range [a, b) or [a, b] depending on rounding." 3574adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return a + (b-a) * self.random() 3584adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3594adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- triangular -------------------- 3604adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3614adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def triangular(self, low=0.0, high=1.0, mode=None): 3624adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Triangular distribution. 3634adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3644adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Continuous distribution bounded by given lower and upper limits, 3654adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao and having a given mode value in-between. 3664adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3674adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao http://en.wikipedia.org/wiki/Triangular_distribution 3684adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3694adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 3704adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao u = self.random() 3714adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao c = 0.5 if mode is None else (mode - low) / (high - low) 3724adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if u > c: 3734adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao u = 1.0 - u 3744adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao c = 1.0 - c 3754adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao low, high = high, low 3764adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return low + (high - low) * (u * c) ** 0.5 3774adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3784adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- normal distribution -------------------- 3794adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3804adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def normalvariate(self, mu, sigma): 3814adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Normal distribution. 3824adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3834adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao mu is the mean, and sigma is the standard deviation. 3844adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3854adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 3864adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # mu = mean, sigma = standard deviation 3874adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3884adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Uses Kinderman and Monahan method. Reference: Kinderman, 3894adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # A.J. and Monahan, J.F., "Computer generation of random 3904adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # variables using the ratio of uniform deviates", ACM Trans 3914adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Math Software, 3, (1977), pp257-260. 3924adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 3934adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao random = self.random 3944adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao while 1: 3954adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao u1 = random() 3964adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao u2 = 1.0 - random() 3974adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao z = NV_MAGICCONST*(u1-0.5)/u2 3984adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao zz = z*z/4.0 3994adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if zz <= -_log(u2): 4004adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao break 4014adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return mu + z*sigma 4024adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4034adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- lognormal distribution -------------------- 4044adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4054adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def lognormvariate(self, mu, sigma): 4064adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Log normal distribution. 4074adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4084adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao If you take the natural logarithm of this distribution, you'll get a 4094adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao normal distribution with mean mu and standard deviation sigma. 4104adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao mu can have any value, and sigma must be greater than zero. 4114adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4124adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 4134adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return _exp(self.normalvariate(mu, sigma)) 4144adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4154adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- exponential distribution -------------------- 4164adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4174adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def expovariate(self, lambd): 4184adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Exponential distribution. 4194adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4204adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao lambd is 1.0 divided by the desired mean. It should be 4214adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao nonzero. (The parameter would be called "lambda", but that is 4224adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao a reserved word in Python.) Returned values range from 0 to 4234adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao positive infinity if lambd is positive, and from negative 4244adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao infinity to 0 if lambd is negative. 4254adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4264adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 4274adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # lambd: rate lambd = 1/mean 4284adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # ('lambda' is a Python reserved word) 4294adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4304adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # we use 1-random() instead of random() to preclude the 4314adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # possibility of taking the log of zero. 4324adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return -_log(1.0 - self.random())/lambd 4334adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4344adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- von Mises distribution -------------------- 4354adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4364adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def vonmisesvariate(self, mu, kappa): 4374adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Circular data distribution. 4384adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4394adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao mu is the mean angle, expressed in radians between 0 and 2*pi, and 4404adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao kappa is the concentration parameter, which must be greater than or 4414adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao equal to zero. If kappa is equal to zero, this distribution reduces 4424adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao to a uniform random angle over the range 0 to 2*pi. 4434adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4444adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 4454adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # mu: mean angle (in radians between 0 and 2*pi) 4464adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # kappa: concentration parameter kappa (>= 0) 4474adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # if kappa = 0 generate uniform random angle 4484adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4494adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Based upon an algorithm published in: Fisher, N.I., 4504adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # "Statistical Analysis of Circular Data", Cambridge 4514adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # University Press, 1993. 4524adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4534adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Thanks to Magnus Kessler for a correction to the 4544adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # implementation of step 4. 4554adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4564adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao random = self.random 4574adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if kappa <= 1e-6: 4584adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return TWOPI * random() 4594adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4604adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao s = 0.5 / kappa 4614adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao r = s + _sqrt(1.0 + s * s) 4624adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4634adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao while 1: 4644adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao u1 = random() 4654adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao z = _cos(_pi * u1) 4664adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4674adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao d = z / (r + z) 4684adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao u2 = random() 4694adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if u2 < 1.0 - d * d or u2 <= (1.0 - d) * _exp(d): 4704adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao break 4714adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4724adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao q = 1.0 / r 4734adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao f = (q + z) / (1.0 + q * z) 4744adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao u3 = random() 4754adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if u3 > 0.5: 4764adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao theta = (mu + _acos(f)) % TWOPI 4774adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao else: 4784adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao theta = (mu - _acos(f)) % TWOPI 4794adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4804adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return theta 4814adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4824adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- gamma distribution -------------------- 4834adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4844adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def gammavariate(self, alpha, beta): 4854adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Gamma distribution. Not the gamma function! 4864adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4874adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Conditions on the parameters are alpha > 0 and beta > 0. 4884adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4894adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao The probability distribution function is: 4904adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4914adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao x ** (alpha - 1) * math.exp(-x / beta) 4924adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao pdf(x) = -------------------------------------- 4934adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao math.gamma(alpha) * beta ** alpha 4944adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4954adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 4964adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4974adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 4984adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 4994adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Warning: a few older sources define the gamma distribution in terms 5004adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # of alpha > -1.0 5014adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if alpha <= 0.0 or beta <= 0.0: 5024adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise ValueError, 'gammavariate: alpha and beta must be > 0.0' 5034adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5044adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao random = self.random 5054adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if alpha > 1.0: 5064adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5074adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Uses R.C.H. Cheng, "The generation of Gamma 5084adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # variables with non-integral shape parameters", 5094adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Applied Statistics, (1977), 26, No. 1, p71-74 5104adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5114adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao ainv = _sqrt(2.0 * alpha - 1.0) 5124adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao bbb = alpha - LOG4 5134adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao ccc = alpha + ainv 5144adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5154adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao while 1: 5164adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao u1 = random() 5174adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if not 1e-7 < u1 < .9999999: 5184adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao continue 5194adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao u2 = 1.0 - random() 5204adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao v = _log(u1/(1.0-u1))/ainv 5214adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao x = alpha*_exp(v) 5224adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao z = u1*u1*u2 5234adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao r = bbb+ccc*v-x 5244adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 5254adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return x * beta 5264adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5274adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao elif alpha == 1.0: 5284adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # expovariate(1) 5294adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao u = random() 5304adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao while u <= 1e-7: 5314adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao u = random() 5324adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return -_log(u) * beta 5334adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5344adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao else: # alpha is between 0 and 1 (exclusive) 5354adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5364adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 5374adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5384adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao while 1: 5394adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao u = random() 5404adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao b = (_e + alpha)/_e 5414adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao p = b*u 5424adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if p <= 1.0: 5434adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao x = p ** (1.0/alpha) 5444adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao else: 5454adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao x = -_log((b-p)/alpha) 5464adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao u1 = random() 5474adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if p > 1.0: 5484adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if u1 <= x ** (alpha - 1.0): 5494adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao break 5504adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao elif u1 <= _exp(-x): 5514adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao break 5524adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return x * beta 5534adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5544adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- Gauss (faster alternative) -------------------- 5554adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5564adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def gauss(self, mu, sigma): 5574adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Gaussian distribution. 5584adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5594adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao mu is the mean, and sigma is the standard deviation. This is 5604adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao slightly faster than the normalvariate() function. 5614adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5624adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Not thread-safe without a lock around calls. 5634adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5644adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 5654adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5664adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # When x and y are two variables from [0, 1), uniformly 5674adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # distributed, then 5684adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # 5694adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # cos(2*pi*x)*sqrt(-2*log(1-y)) 5704adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # sin(2*pi*x)*sqrt(-2*log(1-y)) 5714adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # 5724adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # are two *independent* variables with normal distribution 5734adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # (mu = 0, sigma = 1). 5744adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # (Lambert Meertens) 5754adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # (corrected version; bug discovered by Mike Miller, fixed by LM) 5764adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5774adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Multithreading note: When two threads call this function 5784adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # simultaneously, it is possible that they will receive the 5794adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # same return value. The window is very small though. To 5804adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # avoid this, you have to use a lock around all calls. (I 5814adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # didn't want to slow this down in the serial case by using a 5824adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # lock here.) 5834adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5844adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao random = self.random 5854adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao z = self.gauss_next 5864adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao self.gauss_next = None 5874adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if z is None: 5884adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao x2pi = random() * TWOPI 5894adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao g2rad = _sqrt(-2.0 * _log(1.0 - random())) 5904adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao z = _cos(x2pi) * g2rad 5914adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao self.gauss_next = _sin(x2pi) * g2rad 5924adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5934adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return mu + z*sigma 5944adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 5954adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- beta -------------------- 5964adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## See 5974adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## http://mail.python.org/pipermail/python-bugs-list/2001-January/003752.html 5984adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## for Ivan Frohne's insightful analysis of why the original implementation: 5994adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## 6004adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## def betavariate(self, alpha, beta): 6014adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## # Discrete Event Simulation in C, pp 87-88. 6024adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## 6034adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## y = self.expovariate(alpha) 6044adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## z = self.expovariate(1.0/beta) 6054adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## return z/(y+z) 6064adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## 6074adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## was dead wrong, and how it probably got that way. 6084adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6094adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def betavariate(self, alpha, beta): 6104adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Beta distribution. 6114adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6124adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Conditions on the parameters are alpha > 0 and beta > 0. 6134adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Returned values range between 0 and 1. 6144adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6154adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 6164adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6174adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # This version due to Janne Sinkkonen, and matches all the std 6184adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 6194adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao y = self.gammavariate(alpha, 1.) 6204adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if y == 0: 6214adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return 0.0 6224adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao else: 6234adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return y / (y + self.gammavariate(beta, 1.)) 6244adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6254adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- Pareto -------------------- 6264adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6274adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def paretovariate(self, alpha): 6284adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Pareto distribution. alpha is the shape parameter.""" 6294adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Jain, pg. 495 6304adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6314adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao u = 1.0 - self.random() 6324adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return 1.0 / pow(u, 1.0/alpha) 6334adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6344adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- Weibull -------------------- 6354adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6364adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def weibullvariate(self, alpha, beta): 6374adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Weibull distribution. 6384adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6394adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao alpha is the scale parameter and beta is the shape parameter. 6404adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6414adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 6424adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Jain, pg. 499; bug fix courtesy Bill Arms 6434adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6444adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao u = 1.0 - self.random() 6454adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return alpha * pow(-_log(u), 1.0/beta) 6464adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6474adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- Wichmann-Hill ------------------- 6484adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6494adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaoclass WichmannHill(Random): 6504adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6514adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao VERSION = 1 # used by getstate/setstate 6524adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6534adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def seed(self, a=None): 6544adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Initialize internal state from hashable object. 6554adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6564adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao None or no argument seeds from current time or from an operating 6574adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao system specific randomness source if available. 6584adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6594adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao If a is not None or an int or long, hash(a) is used instead. 6604adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6614adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao If a is an int or long, a is used directly. Distinct values between 6624adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 0 and 27814431486575L inclusive are guaranteed to yield distinct 6634adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao internal states (this guarantee is specific to the default 6644adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Wichmann-Hill generator). 6654adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 6664adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6674adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if a is None: 6684adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao try: 6694adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao a = long(_hexlify(_urandom(16)), 16) 6704adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao except NotImplementedError: 6714adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao import time 6724adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao a = long(time.time() * 256) # use fractional seconds 6734adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6744adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if not isinstance(a, (int, long)): 6754adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao a = hash(a) 6764adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6774adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao a, x = divmod(a, 30268) 6784adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao a, y = divmod(a, 30306) 6794adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao a, z = divmod(a, 30322) 6804adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao self._seed = int(x)+1, int(y)+1, int(z)+1 6814adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6824adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao self.gauss_next = None 6834adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6844adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def random(self): 6854adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Get the next random number in the range [0.0, 1.0).""" 6864adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 6874adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Wichman-Hill random number generator. 6884adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # 6894adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Wichmann, B. A. & Hill, I. D. (1982) 6904adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Algorithm AS 183: 6914adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # An efficient and portable pseudo-random number generator 6924adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Applied Statistics 31 (1982) 188-190 6934adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # 6944adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # see also: 6954adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Correction to Algorithm AS 183 6964adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Applied Statistics 33 (1984) 123 6974adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # 6984adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # McLeod, A. I. (1985) 6994adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # A remark on Algorithm AS 183 7004adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Applied Statistics 34 (1985),198-200 7014adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7024adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # This part is thread-unsafe: 7034adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # BEGIN CRITICAL SECTION 7044adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao x, y, z = self._seed 7054adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao x = (171 * x) % 30269 7064adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao y = (172 * y) % 30307 7074adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao z = (170 * z) % 30323 7084adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao self._seed = x, y, z 7094adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # END CRITICAL SECTION 7104adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7114adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Note: on a platform using IEEE-754 double arithmetic, this can 7124adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # never return 0.0 (asserted by Tim; proof too long for a comment). 7134adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0 7144adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7154adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def getstate(self): 7164adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Return internal state; can be passed to setstate() later.""" 7174adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return self.VERSION, self._seed, self.gauss_next 7184adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7194adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def setstate(self, state): 7204adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Restore internal state from object returned by getstate().""" 7214adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao version = state[0] 7224adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if version == 1: 7234adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao version, self._seed, self.gauss_next = state 7244adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao else: 7254adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise ValueError("state with version %s passed to " 7264adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao "Random.setstate() of version %s" % 7274adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao (version, self.VERSION)) 7284adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7294adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def jumpahead(self, n): 7304adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Act as if n calls to random() were made, but quickly. 7314adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7324adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao n is an int, greater than or equal to 0. 7334adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7344adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Example use: If you have 2 threads and know that each will 7354adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao consume no more than a million random numbers, create two Random 7364adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao objects r1 and r2, then do 7374adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao r2.setstate(r1.getstate()) 7384adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao r2.jumpahead(1000000) 7394adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Then r1 and r2 will use guaranteed-disjoint segments of the full 7404adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao period. 7414adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 7424adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7434adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if not n >= 0: 7444adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise ValueError("n must be >= 0") 7454adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao x, y, z = self._seed 7464adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao x = int(x * pow(171, n, 30269)) % 30269 7474adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao y = int(y * pow(172, n, 30307)) % 30307 7484adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao z = int(z * pow(170, n, 30323)) % 30323 7494adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao self._seed = x, y, z 7504adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7514adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def __whseed(self, x=0, y=0, z=0): 7524adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Set the Wichmann-Hill seed from (x, y, z). 7534adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7544adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao These must be integers in the range [0, 256). 7554adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 7564adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7574adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if not type(x) == type(y) == type(z) == int: 7584adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise TypeError('seeds must be integers') 7594adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): 7604adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise ValueError('seeds must be in range(0, 256)') 7614adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if 0 == x == y == z: 7624adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Initialize from current time 7634adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao import time 7644adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao t = long(time.time() * 256) 7654adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao t = int((t&0xffffff) ^ (t>>24)) 7664adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao t, x = divmod(t, 256) 7674adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao t, y = divmod(t, 256) 7684adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao t, z = divmod(t, 256) 7694adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao # Zero is a poor seed, so substitute 1 7704adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao self._seed = (x or 1, y or 1, z or 1) 7714adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7724adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao self.gauss_next = None 7734adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7744adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def whseed(self, a=None): 7754adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Seed from hashable object's hash code. 7764adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7774adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao None or no argument seeds from current time. It is not guaranteed 7784adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao that objects with distinct hash codes lead to distinct internal 7794adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao states. 7804adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7814adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao This is obsolete, provided for compatibility with the seed routine 7824adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao used prior to Python 2.1. Use the .seed() method instead. 7834adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 7844adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7854adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if a is None: 7864adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao self.__whseed() 7874adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return 7884adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao a = hash(a) 7894adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao a, x = divmod(a, 256) 7904adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao a, y = divmod(a, 256) 7914adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao a, z = divmod(a, 256) 7924adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao x = (x + a) % 256 or 1 7934adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao y = (y + a) % 256 or 1 7944adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao z = (z + a) % 256 or 1 7954adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao self.__whseed(x, y, z) 7964adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7974adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## --------------- Operating System Random Source ------------------ 7984adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 7994adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaoclass SystemRandom(Random): 8004adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Alternate random number generator using sources provided 8014adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao by the operating system (such as /dev/urandom on Unix or 8024adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao CryptGenRandom on Windows). 8034adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 8044adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao Not available on all systems (see os.urandom() for details). 8054adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """ 8064adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 8074adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def random(self): 8084adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """Get the next random number in the range [0.0, 1.0).""" 8094adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return (long(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF 8104adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 8114adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def getrandbits(self, k): 8124adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao """getrandbits(k) -> x. Generates a long int with k random bits.""" 8134adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if k <= 0: 8144adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise ValueError('number of bits must be greater than zero') 8154adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao if k != int(k): 8164adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise TypeError('number of bits should be an integer') 8174adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao bytes = (k + 7) // 8 # bits / 8 and rounded up 8184adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao x = long(_hexlify(_urandom(bytes)), 16) 8194adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return x >> (bytes * 8 - k) # trim excess bits 8204adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 8214adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def _stub(self, *args, **kwds): 8224adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao "Stub method. Not used for a system random number generator." 8234adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao return None 8244adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao seed = jumpahead = _stub 8254adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 8264adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao def _notimplemented(self, *args, **kwds): 8274adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao "Method should not be called for a system random number generator." 8284adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao raise NotImplementedError('System entropy source does not have state.') 8294adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao getstate = setstate = _notimplemented 8304adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 8314adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao## -------------------- test program -------------------- 8324adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 8334adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaodef _test_generator(n, func, args): 8344adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao import time 8354adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao print n, 'times', func.__name__ 8364adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao total = 0.0 8374adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao sqsum = 0.0 8384adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao smallest = 1e10 8394adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao largest = -1e10 8404adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao t0 = time.time() 8414adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao for i in range(n): 8424adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao x = func(*args) 8434adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao total += x 8444adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao sqsum = sqsum + x*x 8454adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao smallest = min(x, smallest) 8464adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao largest = max(x, largest) 8474adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao t1 = time.time() 8484adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao print round(t1-t0, 3), 'sec,', 8494adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao avg = total/n 8504adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao stddev = _sqrt(sqsum/n - avg*avg) 8514adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao print 'avg %g, stddev %g, min %g, max %g' % \ 8524adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao (avg, stddev, smallest, largest) 8534adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 8544adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 8554adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaodef _test(N=2000): 8564adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test_generator(N, random, ()) 8574adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test_generator(N, normalvariate, (0.0, 1.0)) 8584adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test_generator(N, lognormvariate, (0.0, 1.0)) 8594adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test_generator(N, vonmisesvariate, (0.0, 1.0)) 8604adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test_generator(N, gammavariate, (0.01, 1.0)) 8614adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test_generator(N, gammavariate, (0.1, 1.0)) 8624adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test_generator(N, gammavariate, (0.1, 2.0)) 8634adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test_generator(N, gammavariate, (0.5, 1.0)) 8644adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test_generator(N, gammavariate, (0.9, 1.0)) 8654adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test_generator(N, gammavariate, (1.0, 1.0)) 8664adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test_generator(N, gammavariate, (2.0, 1.0)) 8674adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test_generator(N, gammavariate, (20.0, 1.0)) 8684adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test_generator(N, gammavariate, (200.0, 1.0)) 8694adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test_generator(N, gauss, (0.0, 1.0)) 8704adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test_generator(N, betavariate, (3.0, 3.0)) 8714adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test_generator(N, triangular, (0.0, 1.0, 1.0/3.0)) 8724adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 8734adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao# Create one instance, seeded from current time, and export its methods 8744adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao# as module-level functions. The functions share state across all uses 8754adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao#(both in the user's code and in the Python libraries), but that's fine 8764adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao# for most programs and is easier for the casual user than making them 8774adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao# instantiate their own Random() instance. 8784adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 8794adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao_inst = Random() 8804adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaoseed = _inst.seed 8814adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaorandom = _inst.random 8824adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaouniform = _inst.uniform 8834adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaotriangular = _inst.triangular 8844adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaorandint = _inst.randint 8854adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaochoice = _inst.choice 8864adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaorandrange = _inst.randrange 8874adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaosample = _inst.sample 8884adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaoshuffle = _inst.shuffle 8894adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaonormalvariate = _inst.normalvariate 8904adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaolognormvariate = _inst.lognormvariate 8914adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaoexpovariate = _inst.expovariate 8924adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaovonmisesvariate = _inst.vonmisesvariate 8934adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaogammavariate = _inst.gammavariate 8944adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaogauss = _inst.gauss 8954adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaobetavariate = _inst.betavariate 8964adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaoparetovariate = _inst.paretovariate 8974adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaoweibullvariate = _inst.weibullvariate 8984adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaogetstate = _inst.getstate 8994adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaosetstate = _inst.setstate 9004adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaojumpahead = _inst.jumpahead 9014adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaogetrandbits = _inst.getrandbits 9024adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao 9034adfde8bc82dd39f59e0445588c3e599ada477dJosh Gaoif __name__ == '__main__': 9044adfde8bc82dd39f59e0445588c3e599ada477dJosh Gao _test() 905