random.py revision 85e2e4742d0a1accecd02058a7907df36308297e
1e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum"""Random variable generators. 2e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum 3d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters integers 4d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters -------- 5d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters uniform within range 6d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 7d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters sequences 8d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters --------- 9d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters pick random element 10d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters generate random permutation 11d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 12e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum distributions on the real line: 13e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum ------------------------------ 14d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters uniform 15e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum normal (Gaussian) 16e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum lognormal 17e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum negative exponential 18e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum gamma 19e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum beta 20e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum 21e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum distributions on the circle (angles 0 to 2pi) 22e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum --------------------------------------------- 23e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum circular uniform 24e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum von Mises 25e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum 26e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van RossumTranslated from anonymously contributed C/C++ source. 27e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum 28e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van RossumMulti-threading note: the random number generator used here is not 29e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossumthread-safe; it is possible that two calls return the same random 30cd804108548e1939bd8646634ed52ef388ee9f44Tim Petersvalue. But you can instantiate a different instance of Random() in 31cd804108548e1939bd8646634ed52ef388ee9f44Tim Peterseach thread to get generators that don't share state, then use 32cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters.setstate() and .jumpahead() to move the generators to disjoint 33cd804108548e1939bd8646634ed52ef388ee9f44Tim Peterssegments of the full period. 34e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum""" 35d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters# XXX The docstring sucks. 36d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum 37d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersfrom math import log as _log, exp as _exp, pi as _pi, e as _e 38d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersfrom math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin 39d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 40d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersdef _verify(name, expected): 41d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters computed = eval(name) 42d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if abs(computed - expected) > 1e-7: 43d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError( 44d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters "computed value for %s deviates too much " 45d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters "(computed %g, expected %g)" % (name, computed, expected)) 46ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 47d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersNV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) 48d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_verify('NV_MAGICCONST', 1.71552776992141) 4933d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 50d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersTWOPI = 2.0*_pi 51d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_verify('TWOPI', 6.28318530718) 520c9886d589ddebf32de0ca3f027a173222ed383aTim Peters 53d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersLOG4 = _log(4.0) 54d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_verify('LOG4', 1.38629436111989) 550c9886d589ddebf32de0ca3f027a173222ed383aTim Peters 56d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersSG_MAGICCONST = 1.0 + _log(4.5) 57d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_verify('SG_MAGICCONST', 2.50407739677627) 5833d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 59d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersdel _verify 6033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 61d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters# Translated by Guido van Rossum from C source provided by 62d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters# Adrian Baddeley. 6333d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 64d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersclass Random: 6533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 66d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters VERSION = 1 # used by getstate/setstate 6733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 68d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def __init__(self, x=None): 69d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Initialize an instance. 7033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 71d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional argument x controls seeding, as for Random.seed(). 72d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 7333d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 74d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.seed(x) 75d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = None 7633d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 77cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- core generator ------------------- 78cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 79d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Specific to Wichmann-Hill generator. Subclasses wishing to use a 80d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters # different core generator should override the seed(), random(), 81cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # getstate(), setstate() and jumpahead() methods. 82ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 83d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def __whseed(self, x=0, y=0, z=0): 84d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Set the Wichmann-Hill seed from (x, y, z). 85d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 86d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters These must be integers in the range [0, 256). 87d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 88d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 89d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not type(x) == type(y) == type(z) == type(0): 90d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise TypeError('seeds must be integers') 91d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): 92d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError('seeds must be in range(0, 256)') 93d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if 0 == x == y == z: 94d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Initialize from current time 95d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters import time 96d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters t = long(time.time()) * 256 97d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters t = int((t&0xffffff) ^ (t>>24)) 98d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters t, x = divmod(t, 256) 99d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters t, y = divmod(t, 256) 100d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters t, z = divmod(t, 256) 101d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Zero is a poor seed, so substitute 1 102d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self._seed = (x or 1, y or 1, z or 1) 103d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 104cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def random(self): 105cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters """Get the next random number in the range [0.0, 1.0).""" 106cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 107cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Wichman-Hill random number generator. 108cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # 109cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Wichmann, B. A. & Hill, I. D. (1982) 110cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Algorithm AS 183: 111cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # An efficient and portable pseudo-random number generator 112cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Applied Statistics 31 (1982) 188-190 113cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # 114cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # see also: 115cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Correction to Algorithm AS 183 116cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Applied Statistics 33 (1984) 123 117cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # 118cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # McLeod, A. I. (1985) 119cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # A remark on Algorithm AS 183 120cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Applied Statistics 34 (1985),198-200 121cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 122cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # This part is thread-unsafe: 123cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # BEGIN CRITICAL SECTION 124cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters x, y, z = self._seed 125cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters x = (171 * x) % 30269 126cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters y = (172 * y) % 30307 127cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters z = (170 * z) % 30323 128cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters self._seed = x, y, z 129cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # END CRITICAL SECTION 130cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 131cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Note: on a platform using IEEE-754 double arithmetic, this can 132cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # never return 0.0 (asserted by Tim; proof too long for a comment). 133cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0 134cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 135d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def seed(self, a=None): 136cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters """Seed from hashable object's hash code. 137d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 138cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters None or no argument seeds from current time. It is not guaranteed 139cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters that objects with distinct hash codes lead to distinct internal 140cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters states. 141d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 142d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 143d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if a is None: 144d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.__whseed() 145d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters return 146d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters a = hash(a) 147d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters a, x = divmod(a, 256) 148d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters a, y = divmod(a, 256) 149d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters a, z = divmod(a, 256) 150d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = (x + a) % 256 or 1 151d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters y = (y + a) % 256 or 1 152d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = (z + a) % 256 or 1 153d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.__whseed(x, y, z) 154d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 155d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def getstate(self): 156d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Return internal state; can be passed to setstate() later.""" 157d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return self.VERSION, self._seed, self.gauss_next 158d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 159d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def setstate(self, state): 160d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Restore internal state from object returned by getstate().""" 161d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters version = state[0] 162d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if version == 1: 163d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters version, self._seed, self.gauss_next = state 164d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 165d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError("state with version %s passed to " 166d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters "Random.setstate() of version %s" % 167d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters (version, self.VERSION)) 168d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 169d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters def jumpahead(self, n): 170d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters """Act as if n calls to random() were made, but quickly. 171d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters 172d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters n is an int, greater than or equal to 0. 173d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters 174d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters Example use: If you have 2 threads and know that each will 175d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters consume no more than a million random numbers, create two Random 176d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters objects r1 and r2, then do 177d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters r2.setstate(r1.getstate()) 178d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters r2.jumpahead(1000000) 179d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters Then r1 and r2 will use guaranteed-disjoint segments of the full 180d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters period. 181d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters """ 182d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters 183d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters if not n >= 0: 184d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters raise ValueError("n must be >= 0") 185d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters x, y, z = self._seed 186d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters x = int(x * pow(171, n, 30269)) % 30269 187d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters y = int(y * pow(172, n, 30307)) % 30307 188d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters z = int(z * pow(170, n, 30323)) % 30323 189d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters self._seed = x, y, z 190d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters 191cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- Methods below this point do not need to be overridden when 192cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- subclassing for the purpose of using a different core generator. 193d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 194cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- pickle support ------------------- 195d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 196cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __getstate__(self): # for pickle 197cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters return self.getstate() 198d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 199cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __setstate__(self, state): # for pickle 200cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters self.setstate(state) 201cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 202cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- integer methods ------------------- 203d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 204d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def randrange(self, start, stop=None, step=1, int=int, default=None): 205d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random item from range(start, stop[, step]). 206d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 207d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters This fixes the problem with randint() which includes the 208d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters endpoint; in Python this is usually not what you want. 209d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Do not supply the 'int' and 'default' arguments. 210d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 211d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 212d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # This code is a bit messy to make it fast for the 213d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # common case while still doing adequate error checking 214d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istart = int(start) 215d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart != start: 216d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer arg 1 for randrange()" 217d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if stop is default: 218d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart > 0: 219d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return int(self.random() * istart) 220d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 221d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istop = int(stop) 222d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istop != stop: 223d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer stop for randrange()" 224d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if step == 1: 225d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart < istop: 226d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return istart + int(self.random() * 227d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters (istop - istart)) 228d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 229d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istep = int(step) 230d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep != step: 231d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer step for randrange()" 232d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep > 0: 233d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters n = (istop - istart + istep - 1) / istep 234d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif istep < 0: 235d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters n = (istop - istart + istep + 1) / istep 236d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 237d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "zero step for randrange()" 238d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 239d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if n <= 0: 240d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 241d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return istart + istep*int(self.random() * n) 242d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 243d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def randint(self, a, b): 244cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters """Return random integer in range [a, b], including both end points. 245d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 246cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters (Deprecated; use randrange(a, b+1).) 247d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 248d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 249d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return self.randrange(a, b+1) 250d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 251cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- sequence methods ------------------- 252cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 253d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def choice(self, seq): 254d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random element from a non-empty sequence.""" 255d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return seq[int(self.random() * len(seq))] 256d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 257d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def shuffle(self, x, random=None, int=int): 258d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """x, random=random.random -> shuffle list x in place; return None. 259d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 260d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional arg random is a 0-argument function returning a random 261d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters float in [0.0, 1.0); by default, the standard random.random. 262d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 263d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Note that for even rather small len(x), the total number of 264d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters permutations of x is larger than the period of most random number 265d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters generators; this implies that "most" permutations of a long 266d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters sequence can never be generated. 267d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 268d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 269d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if random is None: 270d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 271d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters for i in xrange(len(x)-1, 0, -1): 272cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # pick an element in x[:i+1] with which to exchange x[i] 273d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters j = int(random() * (i+1)) 274d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x[i], x[j] = x[j], x[i] 275d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 276cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- real-valued distributions ------------------- 277cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 278cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- uniform distribution ------------------- 279d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 280d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def uniform(self, a, b): 281d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Get a random number in the range [a, b).""" 282d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return a + (b-a) * self.random() 283ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 284cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- normal distribution -------------------- 285ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 286d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def normalvariate(self, mu, sigma): 287d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mu = mean, sigma = standard deviation 288d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 289d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses Kinderman and Monahan method. Reference: Kinderman, 290d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # A.J. and Monahan, J.F., "Computer generation of random 291d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables using the ratio of uniform deviates", ACM Trans 292d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Math Software, 3, (1977), pp257-260. 293d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 294d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 295d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while 1: 296d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 297d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u2 = random() 298d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = NV_MAGICCONST*(u1-0.5)/u2 299d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters zz = z*z/4.0 300d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if zz <= -_log(u2): 301d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 302d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 303ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 304cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- lognormal distribution -------------------- 305ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 306d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def lognormvariate(self, mu, sigma): 307d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return _exp(self.normalvariate(mu, sigma)) 308ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 309cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- circular uniform -------------------- 310ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 311d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def cunifvariate(self, mean, arc): 312d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mean: mean angle (in radians between 0 and pi) 313d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # arc: range of distribution (in radians between 0 and pi) 314ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 315d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return (mean + arc * (self.random() - 0.5)) % _pi 316ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 317cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- exponential distribution -------------------- 318ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 319d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def expovariate(self, lambd): 320d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lambd: rate lambd = 1/mean 321d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # ('lambda' is a Python reserved word) 322ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 323d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 3240c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 325d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 326d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 327d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return -_log(u)/lambd 328ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 329cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- von Mises distribution -------------------- 330ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 331d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def vonmisesvariate(self, mu, kappa): 332d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mu: mean angle (in radians between 0 and 2*pi) 333d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # kappa: concentration parameter kappa (>= 0) 334d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # if kappa = 0 generate uniform random angle 3355810297052003f28788f6790ac799fe8e5373494Guido van Rossum 336d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Based upon an algorithm published in: Fisher, N.I., 337d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # "Statistical Analysis of Circular Data", Cambridge 338d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # University Press, 1993. 3395810297052003f28788f6790ac799fe8e5373494Guido van Rossum 340d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Thanks to Magnus Kessler for a correction to the 341d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # implementation of step 4. 3425810297052003f28788f6790ac799fe8e5373494Guido van Rossum 343d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 344d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if kappa <= 1e-6: 345d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return TWOPI * random() 346ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 347d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa) 348d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (a - _sqrt(2.0 * a))/(2.0 * kappa) 349d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = (1.0 + b * b)/(2.0 * b) 350ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 351d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while 1: 352d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 353ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 354d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(_pi * u1) 355d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters f = (1.0 + r * z)/(r + z) 356d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters c = kappa * (r - f) 357ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 358d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u2 = random() 359ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 360d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)): 361d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 362ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 363d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u3 = random() 364d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if u3 > 0.5: 365d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters theta = (mu % TWOPI) + _acos(f) 366d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 367d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters theta = (mu % TWOPI) - _acos(f) 368ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 369d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return theta 370ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 371cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- gamma distribution -------------------- 372ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 373d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gammavariate(self, alpha, beta): 374d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # beta times standard gamma 375d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters ainv = _sqrt(2.0 * alpha - 1.0) 376d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return beta * self.stdgamma(alpha, ainv, alpha - LOG4, alpha + ainv) 377d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 378d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def stdgamma(self, alpha, ainv, bbb, ccc): 379d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # ainv = sqrt(2 * alpha - 1) 380d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # bbb = alpha - log(4) 381d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # ccc = alpha + ainv 382d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 383d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 384d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if alpha <= 0.0: 385d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, 'stdgamma: alpha must be > 0.0' 386d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 387d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if alpha > 1.0: 388d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 389d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses R.C.H. Cheng, "The generation of Gamma 390d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables with non-integral shape parameters", 391d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Applied Statistics, (1977), 26, No. 1, p71-74 392d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 393d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while 1: 394d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 395d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u2 = random() 396d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters v = _log(u1/(1.0-u1))/ainv 397d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = alpha*_exp(v) 398d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = u1*u1*u2 399d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = bbb+ccc*v-x 400d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 401d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return x 402d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 403d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif alpha == 1.0: 404d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # expovariate(1) 4050c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 406d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 407d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 408d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return -_log(u) 409d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 410d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: # alpha is between 0 and 1 (exclusive) 411d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 412d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 413d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 414d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while 1: 415d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 416d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (_e + alpha)/_e 417d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters p = b*u 418d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if p <= 1.0: 419d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = pow(p, 1.0/alpha) 420d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 421d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # p > 1 422d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = -_log((b-p)/alpha) 423d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 424d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not (((p <= 1.0) and (u1 > _exp(-x))) or 425d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters ((p > 1) and (u1 > pow(x, alpha - 1.0)))): 426d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 427d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return x 428ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 42995bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 430cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Gauss (faster alternative) -------------------- 43195bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 432d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gauss(self, mu, sigma): 433d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 434d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # When x and y are two variables from [0, 1), uniformly 435d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # distributed, then 436d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 437d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # cos(2*pi*x)*sqrt(-2*log(1-y)) 438d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # sin(2*pi*x)*sqrt(-2*log(1-y)) 439d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 440d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # are two *independent* variables with normal distribution 441d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (mu = 0, sigma = 1). 442d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (Lambert Meertens) 443d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (corrected version; bug discovered by Mike Miller, fixed by LM) 444d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 445d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Multithreading note: When two threads call this function 446d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # simultaneously, it is possible that they will receive the 447d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # same return value. The window is very small though. To 448d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # avoid this, you have to use a lock around all calls. (I 449d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # didn't want to slow this down in the serial case by using a 450d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lock here.) 451d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 452d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 453d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = self.gauss_next 454d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = None 455d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if z is None: 456d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x2pi = random() * TWOPI 457d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters g2rad = _sqrt(-2.0 * _log(1.0 - random())) 458d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(x2pi) * g2rad 459d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = _sin(x2pi) * g2rad 460d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 461d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 46295bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 463cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- beta -------------------- 46485e2e4742d0a1accecd02058a7907df36308297eTim Peters## See 46585e2e4742d0a1accecd02058a7907df36308297eTim Peters## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470 46685e2e4742d0a1accecd02058a7907df36308297eTim Peters## for Ivan Frohne's insightful analysis of why the original implementation: 46785e2e4742d0a1accecd02058a7907df36308297eTim Peters## 46885e2e4742d0a1accecd02058a7907df36308297eTim Peters## def betavariate(self, alpha, beta): 46985e2e4742d0a1accecd02058a7907df36308297eTim Peters## # Discrete Event Simulation in C, pp 87-88. 47085e2e4742d0a1accecd02058a7907df36308297eTim Peters## 47185e2e4742d0a1accecd02058a7907df36308297eTim Peters## y = self.expovariate(alpha) 47285e2e4742d0a1accecd02058a7907df36308297eTim Peters## z = self.expovariate(1.0/beta) 47385e2e4742d0a1accecd02058a7907df36308297eTim Peters## return z/(y+z) 47485e2e4742d0a1accecd02058a7907df36308297eTim Peters## 47585e2e4742d0a1accecd02058a7907df36308297eTim Peters## was dead wrong, and how it probably got that way. 47695bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 477d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def betavariate(self, alpha, beta): 47885e2e4742d0a1accecd02058a7907df36308297eTim Peters # This version due to Janne Sinkkonen, and matches all the std 47985e2e4742d0a1accecd02058a7907df36308297eTim Peters # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 48085e2e4742d0a1accecd02058a7907df36308297eTim Peters y = self.gammavariate(alpha, 1.) 48185e2e4742d0a1accecd02058a7907df36308297eTim Peters if y == 0: 48285e2e4742d0a1accecd02058a7907df36308297eTim Peters return 0.0 48385e2e4742d0a1accecd02058a7907df36308297eTim Peters else: 48485e2e4742d0a1accecd02058a7907df36308297eTim Peters return y / (y + self.gammavariate(beta, 1.)) 48595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 486cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Pareto -------------------- 487cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 488d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def paretovariate(self, alpha): 489d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 495 490cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 491d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = self.random() 492d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return 1.0 / pow(u, 1.0/alpha) 493cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 494cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Weibull -------------------- 495cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 496d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def weibullvariate(self, alpha, beta): 497d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 499; bug fix courtesy Bill Arms 498cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 499d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = self.random() 500d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return alpha * pow(-_log(u), 1.0/beta) 5016c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum 502cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- test program -------------------- 503ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 504d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersdef _test_generator(n, funccall): 5050c9886d589ddebf32de0ca3f027a173222ed383aTim Peters import time 5060c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print n, 'times', funccall 5070c9886d589ddebf32de0ca3f027a173222ed383aTim Peters code = compile(funccall, funccall, 'eval') 5080c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sum = 0.0 5090c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = 0.0 5100c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = 1e10 5110c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = -1e10 5120c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t0 = time.time() 5130c9886d589ddebf32de0ca3f027a173222ed383aTim Peters for i in range(n): 5140c9886d589ddebf32de0ca3f027a173222ed383aTim Peters x = eval(code) 5150c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sum = sum + x 5160c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = sqsum + x*x 5170c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = min(x, smallest) 5180c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = max(x, largest) 5190c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t1 = time.time() 5200c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print round(t1-t0, 3), 'sec,', 5210c9886d589ddebf32de0ca3f027a173222ed383aTim Peters avg = sum/n 522d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters stddev = _sqrt(sqsum/n - avg*avg) 5230c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print 'avg %g, stddev %g, min %g, max %g' % \ 5240c9886d589ddebf32de0ca3f027a173222ed383aTim Peters (avg, stddev, smallest, largest) 525ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 526d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersdef _test(N=200): 527d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters print 'TWOPI =', TWOPI 528d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters print 'LOG4 =', LOG4 529d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters print 'NV_MAGICCONST =', NV_MAGICCONST 530d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters print 'SG_MAGICCONST =', SG_MAGICCONST 531d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'random()') 532d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'normalvariate(0.0, 1.0)') 533d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'lognormvariate(0.0, 1.0)') 534d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'cunifvariate(0.0, 1.0)') 535d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'expovariate(1.0)') 536d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'vonmisesvariate(0.0, 1.0)') 537d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(0.5, 1.0)') 538d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(0.9, 1.0)') 539d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(1.0, 1.0)') 540d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(2.0, 1.0)') 541d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(20.0, 1.0)') 542d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(200.0, 1.0)') 543d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gauss(0.0, 1.0)') 544d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'betavariate(3.0, 3.0)') 545d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'paretovariate(1.0)') 546d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'weibullvariate(1.0, 1.0)') 547d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 548cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Test jumpahead. 549cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters s = getstate() 550cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters jumpahead(N) 551cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters r1 = random() 552cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # now do it the slow way 553cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters setstate(s) 554cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters for i in range(N): 555cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters random() 556cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters r2 = random() 557cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters if r1 != r2: 558cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters raise ValueError("jumpahead test failed " + `(N, r1, r2)`) 559cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 560d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters# Initialize from current time. 561d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_inst = Random() 562d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersseed = _inst.seed 563d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandom = _inst.random 564d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersuniform = _inst.uniform 565d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandint = _inst.randint 566d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterschoice = _inst.choice 567d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandrange = _inst.randrange 568d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersshuffle = _inst.shuffle 569d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersnormalvariate = _inst.normalvariate 570d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterslognormvariate = _inst.lognormvariate 571d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterscunifvariate = _inst.cunifvariate 572d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersexpovariate = _inst.expovariate 573d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersvonmisesvariate = _inst.vonmisesvariate 574d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgammavariate = _inst.gammavariate 575d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersstdgamma = _inst.stdgamma 576d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgauss = _inst.gauss 577d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersbetavariate = _inst.betavariate 578d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersparetovariate = _inst.paretovariate 579d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersweibullvariate = _inst.weibullvariate 580d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgetstate = _inst.getstate 581d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterssetstate = _inst.setstate 582d52269bfd029c4a517ea74c17edd5c3a250c366cTim Petersjumpahead = _inst.jumpahead 583d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 584ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumif __name__ == '__main__': 585d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test() 586