14710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm"""Random variable generators. 24710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 34710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm integers 44710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm -------- 54710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm uniform within range 64710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 74710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm sequences 84710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm --------- 94710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm pick random element 104710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm pick random sample 114710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm generate random permutation 124710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 134710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm distributions on the real line: 144710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm ------------------------------ 154710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm uniform 164710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm triangular 174710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm normal (Gaussian) 184710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm lognormal 194710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm negative exponential 204710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm gamma 214710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm beta 224710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm pareto 234710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Weibull 244710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 254710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm distributions on the circle (angles 0 to 2pi) 264710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm --------------------------------------------- 274710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm circular uniform 284710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm von Mises 294710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 304710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmGeneral notes on the underlying Mersenne Twister core generator: 314710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 324710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm* The period is 2**19937-1. 334710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm* It is one of the most extensively tested generators in existence. 344710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm* Without a direct way to compute N steps forward, the semantics of 354710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm jumpahead(n) are weakened to simply jump to another distant state and rely 364710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm on the large period to avoid overlapping sequences. 374710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm* The random() method is implemented in C, executes in a single Python step, 384710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm and is, therefore, threadsafe. 394710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 404710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm""" 414710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 424710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmfrom __future__ import division 434710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmfrom warnings import warn as _warn 444710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmfrom types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType 454710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmfrom math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil 464710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmfrom math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin 474710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmfrom os import urandom as _urandom 484710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmfrom binascii import hexlify as _hexlify 494710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmimport hashlib as _hashlib 504710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 514710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm__all__ = ["Random","seed","random","uniform","randint","choice","sample", 524710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm "randrange","shuffle","normalvariate","lognormvariate", 534710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm "expovariate","vonmisesvariate","gammavariate","triangular", 544710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm "gauss","betavariate","paretovariate","weibullvariate", 554710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm "getstate","setstate","jumpahead", "WichmannHill", "getrandbits", 564710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm "SystemRandom"] 574710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 584710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmNV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) 594710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmTWOPI = 2.0*_pi 604710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmLOG4 = _log(4.0) 614710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmSG_MAGICCONST = 1.0 + _log(4.5) 624710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmBPF = 53 # Number of bits in a float 634710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmRECIP_BPF = 2**-BPF 644710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 654710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 664710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm# Translated by Guido van Rossum from C source provided by 674710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm# Adrian Baddeley. Adapted by Raymond Hettinger for use with 684710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm# the Mersenne Twister and os.urandom() core generators. 694710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 704710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmimport _random 714710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 724710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmclass Random(_random.Random): 734710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Random number generator base class used by bound module functions. 744710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 754710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Used to instantiate instances of Random to get generators that don't 764710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm share state. Especially useful for multi-threaded programs, creating 774710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm a different instance of Random for each thread, and using the jumpahead() 784710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm method to ensure that the generated sequences seen by each thread don't 794710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm overlap. 804710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 814710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Class Random can also be subclassed if you want to use a different basic 824710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm generator of your own devising: in that case, override the following 834710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm methods: random(), seed(), getstate(), setstate() and jumpahead(). 844710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Optionally, implement a getrandbits() method so that randrange() can cover 854710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm arbitrarily large ranges. 864710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 874710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 884710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 894710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm VERSION = 3 # used by getstate/setstate 904710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 914710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def __init__(self, x=None): 924710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Initialize an instance. 934710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 944710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Optional argument x controls seeding, as for Random.seed(). 954710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 964710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 974710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm self.seed(x) 984710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm self.gauss_next = None 994710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1004710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def seed(self, a=None): 1014710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Initialize internal state from hashable object. 1024710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1034710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm None or no argument seeds from current time or from an operating 1044710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm system specific randomness source if available. 1054710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1064710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm If a is not None or an int or long, hash(a) is used instead. 1074710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 1084710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1094710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if a is None: 1104710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm try: 1114710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm a = long(_hexlify(_urandom(16)), 16) 1124710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm except NotImplementedError: 1134710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm import time 1144710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm a = long(time.time() * 256) # use fractional seconds 1154710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1164710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm super(Random, self).seed(a) 1174710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm self.gauss_next = None 1184710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1194710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def getstate(self): 1204710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Return internal state; can be passed to setstate() later.""" 1214710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return self.VERSION, super(Random, self).getstate(), self.gauss_next 1224710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1234710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def setstate(self, state): 1244710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Restore internal state from object returned by getstate().""" 1254710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm version = state[0] 1264710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if version == 3: 1274710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm version, internalstate, self.gauss_next = state 1284710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm super(Random, self).setstate(internalstate) 1294710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm elif version == 2: 1304710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm version, internalstate, self.gauss_next = state 1314710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # In version 2, the state was saved as signed ints, which causes 1324710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # inconsistencies between 32/64-bit systems. The state is 1334710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # really unsigned 32-bit ints, so we convert negative ints from 1344710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # version 2 to positive longs for version 3. 1354710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm try: 1364710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm internalstate = tuple( long(x) % (2**32) for x in internalstate ) 1374710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm except ValueError, e: 1384710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise TypeError, e 1394710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm super(Random, self).setstate(internalstate) 1404710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm else: 1414710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise ValueError("state with version %s passed to " 1424710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm "Random.setstate() of version %s" % 1434710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm (version, self.VERSION)) 1444710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1454710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def jumpahead(self, n): 1464710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Change the internal state to one that is likely far away 1474710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm from the current state. This method will not be in Py3.x, 1484710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm so it is better to simply reseed. 1494710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 1504710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # The super.jumpahead() method uses shuffling to change state, 1514710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # so it needs a large and "interesting" n to work with. Here, 1524710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # we use hashing to create a large n for the shuffle. 1534710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm s = repr(n) + repr(self.getstate()) 1544710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm n = int(_hashlib.new('sha512', s).hexdigest(), 16) 1554710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm super(Random, self).jumpahead(n) 1564710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1574710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## ---- Methods below this point do not need to be overridden when 1584710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## ---- subclassing for the purpose of using a different core generator. 1594710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1604710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- pickle support ------------------- 1614710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1624710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def __getstate__(self): # for pickle 1634710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return self.getstate() 1644710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1654710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def __setstate__(self, state): # for pickle 1664710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm self.setstate(state) 1674710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1684710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def __reduce__(self): 1694710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return self.__class__, (), self.getstate() 1704710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1714710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- integer methods ------------------- 1724710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1734710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def randrange(self, start, stop=None, step=1, int=int, default=None, 1744710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm maxwidth=1L<<BPF): 1754710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Choose a random item from range(start, stop[, step]). 1764710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1774710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm This fixes the problem with randint() which includes the 1784710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm endpoint; in Python this is usually not what you want. 1794710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Do not supply the 'int', 'default', and 'maxwidth' arguments. 1804710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 1814710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1824710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # This code is a bit messy to make it fast for the 1834710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # common case while still doing adequate error checking. 1844710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm istart = int(start) 1854710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if istart != start: 1864710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise ValueError, "non-integer arg 1 for randrange()" 1874710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if stop is default: 1884710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if istart > 0: 1894710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if istart >= maxwidth: 1904710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return self._randbelow(istart) 1914710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return int(self.random() * istart) 1924710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise ValueError, "empty range for randrange()" 1934710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 1944710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # stop argument supplied. 1954710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm istop = int(stop) 1964710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if istop != stop: 1974710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise ValueError, "non-integer stop for randrange()" 1984710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm width = istop - istart 1994710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if step == 1 and width > 0: 2004710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Note that 2014710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # int(istart + self.random()*width) 2024710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # instead would be incorrect. For example, consider istart 2034710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # = -2 and istop = 0. Then the guts would be in 2044710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # -2.0 to 0.0 exclusive on both ends (ignoring that random() 2054710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # might return 0.0), and because int() truncates toward 0, the 2064710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # final result would be -1 or 0 (instead of -2 or -1). 2074710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # istart + int(self.random()*width) 2084710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # would also be incorrect, for a subtler reason: the RHS 2094710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # can return a long, and then randrange() would also return 2104710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # a long, but we're supposed to return an int (for backward 2114710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # compatibility). 2124710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2134710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if width >= maxwidth: 2144710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return int(istart + self._randbelow(width)) 2154710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return int(istart + int(self.random()*width)) 2164710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if step == 1: 2174710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width) 2184710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2194710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Non-unit step argument supplied. 2204710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm istep = int(step) 2214710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if istep != step: 2224710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise ValueError, "non-integer step for randrange()" 2234710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if istep > 0: 2244710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm n = (width + istep - 1) // istep 2254710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm elif istep < 0: 2264710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm n = (width + istep + 1) // istep 2274710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm else: 2284710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise ValueError, "zero step for randrange()" 2294710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2304710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if n <= 0: 2314710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise ValueError, "empty range for randrange()" 2324710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2334710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if n >= maxwidth: 2344710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return istart + istep*self._randbelow(n) 2354710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return istart + istep*int(self.random() * n) 2364710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2374710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def randint(self, a, b): 2384710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Return random integer in range [a, b], including both end points. 2394710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 2404710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2414710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return self.randrange(a, b+1) 2424710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2434710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def _randbelow(self, n, _log=_log, int=int, _maxwidth=1L<<BPF, 2444710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType): 2454710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Return a random int in the range [0,n) 2464710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2474710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Handles the case where n has more bits than returned 2484710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm by a single call to the underlying generator. 2494710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 2504710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2514710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm try: 2524710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm getrandbits = self.getrandbits 2534710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm except AttributeError: 2544710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm pass 2554710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm else: 2564710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Only call self.getrandbits if the original random() builtin method 2574710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # has not been overridden or if a new getrandbits() was supplied. 2584710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # This assures that the two methods correspond. 2594710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method: 2604710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm k = int(1.00001 + _log(n-1, 2.0)) # 2**k > n-1 > 2**(k-2) 2614710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm r = getrandbits(k) 2624710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm while r >= n: 2634710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm r = getrandbits(k) 2644710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return r 2654710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if n >= _maxwidth: 2664710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _warn("Underlying random() generator does not supply \n" 2674710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm "enough bits to choose from a population range this large") 2684710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return int(self.random() * n) 2694710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2704710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- sequence methods ------------------- 2714710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2724710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def choice(self, seq): 2734710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Choose a random element from a non-empty sequence.""" 2744710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return seq[int(self.random() * len(seq))] # raises IndexError if seq is empty 2754710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2764710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def shuffle(self, x, random=None, int=int): 2774710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """x, random=random.random -> shuffle list x in place; return None. 2784710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2794710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Optional arg random is a 0-argument function returning a random 2804710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm float in [0.0, 1.0); by default, the standard random.random. 2814710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 2824710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2834710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if random is None: 2844710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm random = self.random 2854710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm for i in reversed(xrange(1, len(x))): 2864710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # pick an element in x[:i+1] with which to exchange x[i] 2874710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm j = int(random() * (i+1)) 2884710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm x[i], x[j] = x[j], x[i] 2894710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2904710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def sample(self, population, k): 2914710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Chooses k unique random elements from a population sequence. 2924710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2934710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Returns a new list containing elements from the population while 2944710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm leaving the original population unchanged. The resulting list is 2954710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm in selection order so that all sub-slices will also be valid random 2964710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm samples. This allows raffle winners (the sample) to be partitioned 2974710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm into grand prize and second place winners (the subslices). 2984710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 2994710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Members of the population need not be hashable or unique. If the 3004710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm population contains repeats, then each occurrence is a possible 3014710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm selection in the sample. 3024710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3034710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm To choose a sample in a range of integers, use xrange as an argument. 3044710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm This is especially fast and space efficient for sampling from a 3054710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm large population: sample(xrange(10000000), 60) 3064710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 3074710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3084710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Sampling without replacement entails tracking either potential 3094710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # selections (the pool) in a list or previous selections in a set. 3104710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3114710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # When the number of selections is small compared to the 3124710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # population, then tracking selections is efficient, requiring 3134710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # only a small set and an occasional reselection. For 3144710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # a larger number of selections, the pool tracking method is 3154710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # preferred since the list takes less space than the 3164710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # set and it doesn't suffer from frequent reselections. 3174710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3184710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm n = len(population) 3194710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if not 0 <= k <= n: 3204710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise ValueError("sample larger than population") 3214710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm random = self.random 3224710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _int = int 3234710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm result = [None] * k 3244710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm setsize = 21 # size of a small set minus size of an empty list 3254710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if k > 5: 3264710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets 3274710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if n <= setsize or hasattr(population, "keys"): 3284710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # An n-length list is smaller than a k-length set, or this is a 3294710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # mapping type so the other algorithm wouldn't work. 3304710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm pool = list(population) 3314710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm for i in xrange(k): # invariant: non-selected at [0,n-i) 3324710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm j = _int(random() * (n-i)) 3334710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm result[i] = pool[j] 3344710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm pool[j] = pool[n-i-1] # move non-selected item into vacancy 3354710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm else: 3364710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm try: 3374710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm selected = set() 3384710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm selected_add = selected.add 3394710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm for i in xrange(k): 3404710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm j = _int(random() * n) 3414710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm while j in selected: 3424710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm j = _int(random() * n) 3434710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm selected_add(j) 3444710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm result[i] = population[j] 3454710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm except (TypeError, KeyError): # handle (at least) sets 3464710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if isinstance(population, list): 3474710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise 3484710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return self.sample(tuple(population), k) 3494710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return result 3504710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3514710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- real-valued distributions ------------------- 3524710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3534710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- uniform distribution ------------------- 3544710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3554710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def uniform(self, a, b): 3564710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm "Get a random number in the range [a, b) or [a, b] depending on rounding." 3574710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return a + (b-a) * self.random() 3584710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3594710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- triangular -------------------- 3604710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3614710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def triangular(self, low=0.0, high=1.0, mode=None): 3624710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Triangular distribution. 3634710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3644710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Continuous distribution bounded by given lower and upper limits, 3654710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm and having a given mode value in-between. 3664710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3674710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm http://en.wikipedia.org/wiki/Triangular_distribution 3684710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3694710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 3704710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u = self.random() 3714710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm c = 0.5 if mode is None else (mode - low) / (high - low) 3724710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if u > c: 3734710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u = 1.0 - u 3744710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm c = 1.0 - c 3754710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm low, high = high, low 3764710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return low + (high - low) * (u * c) ** 0.5 3774710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3784710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- normal distribution -------------------- 3794710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3804710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def normalvariate(self, mu, sigma): 3814710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Normal distribution. 3824710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3834710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm mu is the mean, and sigma is the standard deviation. 3844710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3854710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 3864710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # mu = mean, sigma = standard deviation 3874710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3884710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Uses Kinderman and Monahan method. Reference: Kinderman, 3894710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # A.J. and Monahan, J.F., "Computer generation of random 3904710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # variables using the ratio of uniform deviates", ACM Trans 3914710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Math Software, 3, (1977), pp257-260. 3924710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 3934710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm random = self.random 3944710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm while 1: 3954710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u1 = random() 3964710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u2 = 1.0 - random() 3974710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm z = NV_MAGICCONST*(u1-0.5)/u2 3984710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm zz = z*z/4.0 3994710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if zz <= -_log(u2): 4004710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm break 4014710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return mu + z*sigma 4024710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4034710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- lognormal distribution -------------------- 4044710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4054710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def lognormvariate(self, mu, sigma): 4064710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Log normal distribution. 4074710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4084710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm If you take the natural logarithm of this distribution, you'll get a 4094710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm normal distribution with mean mu and standard deviation sigma. 4104710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm mu can have any value, and sigma must be greater than zero. 4114710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4124710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 4134710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return _exp(self.normalvariate(mu, sigma)) 4144710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4154710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- exponential distribution -------------------- 4164710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4174710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def expovariate(self, lambd): 4184710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Exponential distribution. 4194710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4204710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm lambd is 1.0 divided by the desired mean. It should be 4214710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm nonzero. (The parameter would be called "lambda", but that is 4224710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm a reserved word in Python.) Returned values range from 0 to 4234710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm positive infinity if lambd is positive, and from negative 4244710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm infinity to 0 if lambd is negative. 4254710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4264710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 4274710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # lambd: rate lambd = 1/mean 4284710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # ('lambda' is a Python reserved word) 4294710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4304710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm random = self.random 4314710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u = random() 4324710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm while u <= 1e-7: 4334710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u = random() 4344710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return -_log(u)/lambd 4354710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4364710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- von Mises distribution -------------------- 4374710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4384710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def vonmisesvariate(self, mu, kappa): 4394710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Circular data distribution. 4404710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4414710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm mu is the mean angle, expressed in radians between 0 and 2*pi, and 4424710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm kappa is the concentration parameter, which must be greater than or 4434710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm equal to zero. If kappa is equal to zero, this distribution reduces 4444710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm to a uniform random angle over the range 0 to 2*pi. 4454710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4464710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 4474710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # mu: mean angle (in radians between 0 and 2*pi) 4484710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # kappa: concentration parameter kappa (>= 0) 4494710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # if kappa = 0 generate uniform random angle 4504710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4514710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Based upon an algorithm published in: Fisher, N.I., 4524710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # "Statistical Analysis of Circular Data", Cambridge 4534710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # University Press, 1993. 4544710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4554710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Thanks to Magnus Kessler for a correction to the 4564710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # implementation of step 4. 4574710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4584710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm random = self.random 4594710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if kappa <= 1e-6: 4604710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return TWOPI * random() 4614710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4624710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa) 4634710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm b = (a - _sqrt(2.0 * a))/(2.0 * kappa) 4644710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm r = (1.0 + b * b)/(2.0 * b) 4654710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4664710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm while 1: 4674710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u1 = random() 4684710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4694710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm z = _cos(_pi * u1) 4704710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm f = (1.0 + r * z)/(r + z) 4714710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm c = kappa * (r - f) 4724710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4734710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u2 = random() 4744710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4754710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if u2 < c * (2.0 - c) or u2 <= c * _exp(1.0 - c): 4764710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm break 4774710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4784710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u3 = random() 4794710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if u3 > 0.5: 4804710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm theta = (mu % TWOPI) + _acos(f) 4814710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm else: 4824710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm theta = (mu % TWOPI) - _acos(f) 4834710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4844710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return theta 4854710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4864710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- gamma distribution -------------------- 4874710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4884710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def gammavariate(self, alpha, beta): 4894710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Gamma distribution. Not the gamma function! 4904710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4914710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Conditions on the parameters are alpha > 0 and beta > 0. 4924710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4934710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm The probability distribution function is: 4944710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4954710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm x ** (alpha - 1) * math.exp(-x / beta) 4964710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm pdf(x) = -------------------------------------- 4974710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm math.gamma(alpha) * beta ** alpha 4984710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 4994710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 5004710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5014710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 5024710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5034710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Warning: a few older sources define the gamma distribution in terms 5044710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # of alpha > -1.0 5054710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if alpha <= 0.0 or beta <= 0.0: 5064710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise ValueError, 'gammavariate: alpha and beta must be > 0.0' 5074710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5084710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm random = self.random 5094710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if alpha > 1.0: 5104710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5114710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Uses R.C.H. Cheng, "The generation of Gamma 5124710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # variables with non-integral shape parameters", 5134710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Applied Statistics, (1977), 26, No. 1, p71-74 5144710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5154710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm ainv = _sqrt(2.0 * alpha - 1.0) 5164710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm bbb = alpha - LOG4 5174710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm ccc = alpha + ainv 5184710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5194710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm while 1: 5204710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u1 = random() 5214710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if not 1e-7 < u1 < .9999999: 5224710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm continue 5234710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u2 = 1.0 - random() 5244710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm v = _log(u1/(1.0-u1))/ainv 5254710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm x = alpha*_exp(v) 5264710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm z = u1*u1*u2 5274710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm r = bbb+ccc*v-x 5284710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 5294710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return x * beta 5304710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5314710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm elif alpha == 1.0: 5324710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # expovariate(1) 5334710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u = random() 5344710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm while u <= 1e-7: 5354710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u = random() 5364710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return -_log(u) * beta 5374710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5384710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm else: # alpha is between 0 and 1 (exclusive) 5394710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5404710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 5414710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5424710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm while 1: 5434710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u = random() 5444710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm b = (_e + alpha)/_e 5454710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm p = b*u 5464710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if p <= 1.0: 5474710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm x = p ** (1.0/alpha) 5484710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm else: 5494710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm x = -_log((b-p)/alpha) 5504710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u1 = random() 5514710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if p > 1.0: 5524710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if u1 <= x ** (alpha - 1.0): 5534710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm break 5544710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm elif u1 <= _exp(-x): 5554710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm break 5564710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return x * beta 5574710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5584710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- Gauss (faster alternative) -------------------- 5594710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5604710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def gauss(self, mu, sigma): 5614710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Gaussian distribution. 5624710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5634710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm mu is the mean, and sigma is the standard deviation. This is 5644710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm slightly faster than the normalvariate() function. 5654710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5664710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Not thread-safe without a lock around calls. 5674710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5684710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 5694710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5704710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # When x and y are two variables from [0, 1), uniformly 5714710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # distributed, then 5724710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # 5734710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # cos(2*pi*x)*sqrt(-2*log(1-y)) 5744710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # sin(2*pi*x)*sqrt(-2*log(1-y)) 5754710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # 5764710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # are two *independent* variables with normal distribution 5774710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # (mu = 0, sigma = 1). 5784710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # (Lambert Meertens) 5794710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # (corrected version; bug discovered by Mike Miller, fixed by LM) 5804710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5814710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Multithreading note: When two threads call this function 5824710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # simultaneously, it is possible that they will receive the 5834710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # same return value. The window is very small though. To 5844710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # avoid this, you have to use a lock around all calls. (I 5854710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # didn't want to slow this down in the serial case by using a 5864710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # lock here.) 5874710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5884710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm random = self.random 5894710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm z = self.gauss_next 5904710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm self.gauss_next = None 5914710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if z is None: 5924710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm x2pi = random() * TWOPI 5934710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm g2rad = _sqrt(-2.0 * _log(1.0 - random())) 5944710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm z = _cos(x2pi) * g2rad 5954710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm self.gauss_next = _sin(x2pi) * g2rad 5964710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5974710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return mu + z*sigma 5984710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 5994710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- beta -------------------- 6004710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## See 6014710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## http://mail.python.org/pipermail/python-bugs-list/2001-January/003752.html 6024710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## for Ivan Frohne's insightful analysis of why the original implementation: 6034710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## 6044710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## def betavariate(self, alpha, beta): 6054710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## # Discrete Event Simulation in C, pp 87-88. 6064710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## 6074710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## y = self.expovariate(alpha) 6084710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## z = self.expovariate(1.0/beta) 6094710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## return z/(y+z) 6104710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## 6114710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## was dead wrong, and how it probably got that way. 6124710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6134710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def betavariate(self, alpha, beta): 6144710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Beta distribution. 6154710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6164710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Conditions on the parameters are alpha > 0 and beta > 0. 6174710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Returned values range between 0 and 1. 6184710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6194710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 6204710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6214710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # This version due to Janne Sinkkonen, and matches all the std 6224710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 6234710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm y = self.gammavariate(alpha, 1.) 6244710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if y == 0: 6254710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return 0.0 6264710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm else: 6274710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return y / (y + self.gammavariate(beta, 1.)) 6284710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6294710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- Pareto -------------------- 6304710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6314710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def paretovariate(self, alpha): 6324710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Pareto distribution. alpha is the shape parameter.""" 6334710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Jain, pg. 495 6344710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6354710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u = 1.0 - self.random() 6364710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return 1.0 / pow(u, 1.0/alpha) 6374710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6384710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- Weibull -------------------- 6394710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6404710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def weibullvariate(self, alpha, beta): 6414710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Weibull distribution. 6424710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6434710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm alpha is the scale parameter and beta is the shape parameter. 6444710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6454710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 6464710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Jain, pg. 499; bug fix courtesy Bill Arms 6474710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6484710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm u = 1.0 - self.random() 6494710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return alpha * pow(-_log(u), 1.0/beta) 6504710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6514710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- Wichmann-Hill ------------------- 6524710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6534710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmclass WichmannHill(Random): 6544710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6554710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm VERSION = 1 # used by getstate/setstate 6564710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6574710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def seed(self, a=None): 6584710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Initialize internal state from hashable object. 6594710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6604710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm None or no argument seeds from current time or from an operating 6614710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm system specific randomness source if available. 6624710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6634710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm If a is not None or an int or long, hash(a) is used instead. 6644710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6654710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm If a is an int or long, a is used directly. Distinct values between 6664710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 0 and 27814431486575L inclusive are guaranteed to yield distinct 6674710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm internal states (this guarantee is specific to the default 6684710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Wichmann-Hill generator). 6694710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 6704710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6714710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if a is None: 6724710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm try: 6734710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm a = long(_hexlify(_urandom(16)), 16) 6744710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm except NotImplementedError: 6754710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm import time 6764710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm a = long(time.time() * 256) # use fractional seconds 6774710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6784710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if not isinstance(a, (int, long)): 6794710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm a = hash(a) 6804710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6814710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm a, x = divmod(a, 30268) 6824710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm a, y = divmod(a, 30306) 6834710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm a, z = divmod(a, 30322) 6844710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm self._seed = int(x)+1, int(y)+1, int(z)+1 6854710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6864710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm self.gauss_next = None 6874710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6884710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def random(self): 6894710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Get the next random number in the range [0.0, 1.0).""" 6904710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 6914710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Wichman-Hill random number generator. 6924710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # 6934710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Wichmann, B. A. & Hill, I. D. (1982) 6944710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Algorithm AS 183: 6954710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # An efficient and portable pseudo-random number generator 6964710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Applied Statistics 31 (1982) 188-190 6974710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # 6984710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # see also: 6994710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Correction to Algorithm AS 183 7004710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Applied Statistics 33 (1984) 123 7014710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # 7024710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # McLeod, A. I. (1985) 7034710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # A remark on Algorithm AS 183 7044710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Applied Statistics 34 (1985),198-200 7054710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 7064710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # This part is thread-unsafe: 7074710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # BEGIN CRITICAL SECTION 7084710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm x, y, z = self._seed 7094710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm x = (171 * x) % 30269 7104710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm y = (172 * y) % 30307 7114710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm z = (170 * z) % 30323 7124710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm self._seed = x, y, z 7134710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # END CRITICAL SECTION 7144710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 7154710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Note: on a platform using IEEE-754 double arithmetic, this can 7164710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # never return 0.0 (asserted by Tim; proof too long for a comment). 7174710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0 7184710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 7194710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def getstate(self): 7204710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Return internal state; can be passed to setstate() later.""" 7214710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return self.VERSION, self._seed, self.gauss_next 7224710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 7234710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def setstate(self, state): 7244710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Restore internal state from object returned by getstate().""" 7254710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm version = state[0] 7264710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if version == 1: 7274710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm version, self._seed, self.gauss_next = state 7284710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm else: 7294710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise ValueError("state with version %s passed to " 7304710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm "Random.setstate() of version %s" % 7314710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm (version, self.VERSION)) 7324710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 7334710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def jumpahead(self, n): 7344710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Act as if n calls to random() were made, but quickly. 7354710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 7364710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm n is an int, greater than or equal to 0. 7374710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 7384710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Example use: If you have 2 threads and know that each will 7394710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm consume no more than a million random numbers, create two Random 7404710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm objects r1 and r2, then do 7414710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm r2.setstate(r1.getstate()) 7424710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm r2.jumpahead(1000000) 7434710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Then r1 and r2 will use guaranteed-disjoint segments of the full 7444710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm period. 7454710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 7464710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 7474710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if not n >= 0: 7484710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise ValueError("n must be >= 0") 7494710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm x, y, z = self._seed 7504710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm x = int(x * pow(171, n, 30269)) % 30269 7514710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm y = int(y * pow(172, n, 30307)) % 30307 7524710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm z = int(z * pow(170, n, 30323)) % 30323 7534710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm self._seed = x, y, z 7544710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 7554710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def __whseed(self, x=0, y=0, z=0): 7564710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Set the Wichmann-Hill seed from (x, y, z). 7574710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 7584710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm These must be integers in the range [0, 256). 7594710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 7604710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 7614710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if not type(x) == type(y) == type(z) == int: 7624710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise TypeError('seeds must be integers') 7634710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): 7644710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise ValueError('seeds must be in range(0, 256)') 7654710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if 0 == x == y == z: 7664710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Initialize from current time 7674710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm import time 7684710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm t = long(time.time() * 256) 7694710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm t = int((t&0xffffff) ^ (t>>24)) 7704710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm t, x = divmod(t, 256) 7714710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm t, y = divmod(t, 256) 7724710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm t, z = divmod(t, 256) 7734710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm # Zero is a poor seed, so substitute 1 7744710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm self._seed = (x or 1, y or 1, z or 1) 7754710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 7764710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm self.gauss_next = None 7774710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 7784710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def whseed(self, a=None): 7794710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Seed from hashable object's hash code. 7804710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 7814710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm None or no argument seeds from current time. It is not guaranteed 7824710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm that objects with distinct hash codes lead to distinct internal 7834710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm states. 7844710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 7854710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm This is obsolete, provided for compatibility with the seed routine 7864710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm used prior to Python 2.1. Use the .seed() method instead. 7874710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 7884710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 7894710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if a is None: 7904710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm self.__whseed() 7914710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return 7924710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm a = hash(a) 7934710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm a, x = divmod(a, 256) 7944710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm a, y = divmod(a, 256) 7954710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm a, z = divmod(a, 256) 7964710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm x = (x + a) % 256 or 1 7974710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm y = (y + a) % 256 or 1 7984710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm z = (z + a) % 256 or 1 7994710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm self.__whseed(x, y, z) 8004710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 8014710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## --------------- Operating System Random Source ------------------ 8024710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 8034710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmclass SystemRandom(Random): 8044710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Alternate random number generator using sources provided 8054710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm by the operating system (such as /dev/urandom on Unix or 8064710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm CryptGenRandom on Windows). 8074710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 8084710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm Not available on all systems (see os.urandom() for details). 8094710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """ 8104710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 8114710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def random(self): 8124710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """Get the next random number in the range [0.0, 1.0).""" 8134710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return (long(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF 8144710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 8154710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def getrandbits(self, k): 8164710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm """getrandbits(k) -> x. Generates a long int with k random bits.""" 8174710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if k <= 0: 8184710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise ValueError('number of bits must be greater than zero') 8194710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm if k != int(k): 8204710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise TypeError('number of bits should be an integer') 8214710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm bytes = (k + 7) // 8 # bits / 8 and rounded up 8224710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm x = long(_hexlify(_urandom(bytes)), 16) 8234710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return x >> (bytes * 8 - k) # trim excess bits 8244710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 8254710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def _stub(self, *args, **kwds): 8264710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm "Stub method. Not used for a system random number generator." 8274710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm return None 8284710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm seed = jumpahead = _stub 8294710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 8304710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm def _notimplemented(self, *args, **kwds): 8314710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm "Method should not be called for a system random number generator." 8324710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm raise NotImplementedError('System entropy source does not have state.') 8334710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm getstate = setstate = _notimplemented 8344710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 8354710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm## -------------------- test program -------------------- 8364710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 8374710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmdef _test_generator(n, func, args): 8384710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm import time 8394710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm print n, 'times', func.__name__ 8404710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm total = 0.0 8414710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm sqsum = 0.0 8424710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm smallest = 1e10 8434710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm largest = -1e10 8444710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm t0 = time.time() 8454710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm for i in range(n): 8464710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm x = func(*args) 8474710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm total += x 8484710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm sqsum = sqsum + x*x 8494710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm smallest = min(x, smallest) 8504710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm largest = max(x, largest) 8514710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm t1 = time.time() 8524710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm print round(t1-t0, 3), 'sec,', 8534710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm avg = total/n 8544710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm stddev = _sqrt(sqsum/n - avg*avg) 8554710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm print 'avg %g, stddev %g, min %g, max %g' % \ 8564710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm (avg, stddev, smallest, largest) 8574710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 8584710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 8594710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmdef _test(N=2000): 8604710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test_generator(N, random, ()) 8614710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test_generator(N, normalvariate, (0.0, 1.0)) 8624710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test_generator(N, lognormvariate, (0.0, 1.0)) 8634710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test_generator(N, vonmisesvariate, (0.0, 1.0)) 8644710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test_generator(N, gammavariate, (0.01, 1.0)) 8654710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test_generator(N, gammavariate, (0.1, 1.0)) 8664710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test_generator(N, gammavariate, (0.1, 2.0)) 8674710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test_generator(N, gammavariate, (0.5, 1.0)) 8684710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test_generator(N, gammavariate, (0.9, 1.0)) 8694710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test_generator(N, gammavariate, (1.0, 1.0)) 8704710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test_generator(N, gammavariate, (2.0, 1.0)) 8714710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test_generator(N, gammavariate, (20.0, 1.0)) 8724710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test_generator(N, gammavariate, (200.0, 1.0)) 8734710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test_generator(N, gauss, (0.0, 1.0)) 8744710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test_generator(N, betavariate, (3.0, 3.0)) 8754710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test_generator(N, triangular, (0.0, 1.0, 1.0/3.0)) 8764710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 8774710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm# Create one instance, seeded from current time, and export its methods 8784710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm# as module-level functions. The functions share state across all uses 8794710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm#(both in the user's code and in the Python libraries), but that's fine 8804710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm# for most programs and is easier for the casual user than making them 8814710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm# instantiate their own Random() instance. 8824710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 8834710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm_inst = Random() 8844710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmseed = _inst.seed 8854710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmrandom = _inst.random 8864710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmuniform = _inst.uniform 8874710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmtriangular = _inst.triangular 8884710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmrandint = _inst.randint 8894710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmchoice = _inst.choice 8904710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmrandrange = _inst.randrange 8914710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmsample = _inst.sample 8924710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmshuffle = _inst.shuffle 8934710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmnormalvariate = _inst.normalvariate 8944710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmlognormvariate = _inst.lognormvariate 8954710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmexpovariate = _inst.expovariate 8964710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmvonmisesvariate = _inst.vonmisesvariate 8974710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmgammavariate = _inst.gammavariate 8984710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmgauss = _inst.gauss 8994710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmbetavariate = _inst.betavariate 9004710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmparetovariate = _inst.paretovariate 9014710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmweibullvariate = _inst.weibullvariate 9024710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmgetstate = _inst.getstate 9034710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmsetstate = _inst.setstate 9044710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmjumpahead = _inst.jumpahead 9054710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmgetrandbits = _inst.getrandbits 9064710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm 9074710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylmif __name__ == '__main__': 9084710c53dcad1ebf3755f3efb9e80ac24bd72a9b2darylm _test() 909