random.py revision 0de65807e6bdc5254f5a7e99b2f39adeea6b883b
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 28e360d9507a698290484f2d6b2ff7941db3d30045Tim PetersMulti-threading note: the random number generator used here is not thread- 29e360d9507a698290484f2d6b2ff7941db3d30045Tim Peterssafe; it is possible that two calls return the same random value. However, 30e360d9507a698290484f2d6b2ff7941db3d30045Tim Petersyou can instantiate a different instance of Random() in each thread to get 31e360d9507a698290484f2d6b2ff7941db3d30045Tim Petersgenerators that don't share state, then use .setstate() and .jumpahead() to 32e360d9507a698290484f2d6b2ff7941db3d30045Tim Petersmove the generators to disjoint segments of the full period. For example, 33e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters 34e360d9507a698290484f2d6b2ff7941db3d30045Tim Petersdef create_generators(num, delta, firstseed=None): 35e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters ""\"Return list of num distinct generators. 36e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters Each generator has its own unique segment of delta elements from 37e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters Random.random()'s full period. 38e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters Seed the first generator with optional arg firstseed (default is 39e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters None, to seed from current time). 40e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters ""\" 41e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters 42e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters from random import Random 43e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters g = Random(firstseed) 44e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters result = [g] 45e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters for i in range(num - 1): 46e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters laststate = g.getstate() 47e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters g = Random() 48e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters g.setstate(laststate) 49e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters g.jumpahead(delta) 50e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters result.append(g) 51e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters return result 52e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters 53e360d9507a698290484f2d6b2ff7941db3d30045Tim Petersgens = create_generators(10, 1000000) 54e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters 55e360d9507a698290484f2d6b2ff7941db3d30045Tim PetersThat creates 10 distinct generators, which can be passed out to 10 distinct 56e360d9507a698290484f2d6b2ff7941db3d30045Tim Petersthreads. The generators don't share state so can be called safely in 57e360d9507a698290484f2d6b2ff7941db3d30045Tim Petersparallel. So long as no thread calls its g.random() more than a million 58e360d9507a698290484f2d6b2ff7941db3d30045Tim Peterstimes (the second argument to create_generators), the sequences seen by 59e360d9507a698290484f2d6b2ff7941db3d30045Tim Peterseach thread will not overlap. 60e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters 61e360d9507a698290484f2d6b2ff7941db3d30045Tim PetersThe period of the underlying Wichmann-Hill generator is 6,953,607,871,644, 62e360d9507a698290484f2d6b2ff7941db3d30045Tim Petersand that limits how far this technique can be pushed. 63e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters 64e360d9507a698290484f2d6b2ff7941db3d30045Tim PetersJust for fun, note that since we know the period, .jumpahead() can also be 65e360d9507a698290484f2d6b2ff7941db3d30045Tim Petersused to "move backward in time": 66e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters 67e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters>>> g = Random(42) # arbitrary 68e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters>>> g.random() 690de88fc4b108751b86443852b6741680d704168fTim Peters0.25420336316883324 70e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters>>> g.jumpahead(6953607871644L - 1) # move *back* one 71e360d9507a698290484f2d6b2ff7941db3d30045Tim Peters>>> g.random() 720de88fc4b108751b86443852b6741680d704168fTim Peters0.25420336316883324 73e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum""" 74d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters# XXX The docstring sucks. 75d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum 76d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersfrom math import log as _log, exp as _exp, pi as _pi, e as _e 77d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersfrom math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin 78d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 790de65807e6bdc5254f5a7e99b2f39adeea6b883bSkip Montanaro__all__ = ["Random","seed","random","uniform","randint","choice", 800de65807e6bdc5254f5a7e99b2f39adeea6b883bSkip Montanaro "randrange","shuffle","normalvariate","lognormvariate", 810de65807e6bdc5254f5a7e99b2f39adeea6b883bSkip Montanaro "cunifvariate","expovariate","vonmisesvariate","gammavariate", 820de65807e6bdc5254f5a7e99b2f39adeea6b883bSkip Montanaro "stdgamma","gauss","betavariate","paretovariate","weibullvariate", 830de65807e6bdc5254f5a7e99b2f39adeea6b883bSkip Montanaro "getstate","setstate","jumpahead","whseed"] 840de65807e6bdc5254f5a7e99b2f39adeea6b883bSkip Montanaro 85d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersdef _verify(name, expected): 86d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters computed = eval(name) 87d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if abs(computed - expected) > 1e-7: 88d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError( 89d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters "computed value for %s deviates too much " 90d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters "(computed %g, expected %g)" % (name, computed, expected)) 91ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 92d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersNV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) 93d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_verify('NV_MAGICCONST', 1.71552776992141) 9433d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 95d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersTWOPI = 2.0*_pi 96d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_verify('TWOPI', 6.28318530718) 970c9886d589ddebf32de0ca3f027a173222ed383aTim Peters 98d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersLOG4 = _log(4.0) 99d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_verify('LOG4', 1.38629436111989) 1000c9886d589ddebf32de0ca3f027a173222ed383aTim Peters 101d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersSG_MAGICCONST = 1.0 + _log(4.5) 102d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_verify('SG_MAGICCONST', 2.50407739677627) 10333d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 104d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersdel _verify 10533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 106d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters# Translated by Guido van Rossum from C source provided by 107d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters# Adrian Baddeley. 10833d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 109d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersclass Random: 11033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 111d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters VERSION = 1 # used by getstate/setstate 11233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 113d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def __init__(self, x=None): 114d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Initialize an instance. 11533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 116d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional argument x controls seeding, as for Random.seed(). 117d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 11833d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 119d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.seed(x) 120d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = None 12133d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 122cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- core generator ------------------- 123cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 124d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Specific to Wichmann-Hill generator. Subclasses wishing to use a 125d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters # different core generator should override the seed(), random(), 126cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # getstate(), setstate() and jumpahead() methods. 127ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 1280de88fc4b108751b86443852b6741680d704168fTim Peters def seed(self, a=None): 1290de88fc4b108751b86443852b6741680d704168fTim Peters """Initialize internal state from hashable object. 130d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 1310de88fc4b108751b86443852b6741680d704168fTim Peters None or no argument seeds from current time. 1320de88fc4b108751b86443852b6741680d704168fTim Peters 133bcd725fc456faca13f4598f87c0517f917711cdaTim Peters If a is not None or an int or long, hash(a) is used instead. 1340de88fc4b108751b86443852b6741680d704168fTim Peters 1350de88fc4b108751b86443852b6741680d704168fTim Peters If a is an int or long, a is used directly. Distinct values between 1360de88fc4b108751b86443852b6741680d704168fTim Peters 0 and 27814431486575L inclusive are guaranteed to yield distinct 1370de88fc4b108751b86443852b6741680d704168fTim Peters internal states (this guarantee is specific to the default 1380de88fc4b108751b86443852b6741680d704168fTim Peters Wichmann-Hill generator). 139d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 140d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 1410de88fc4b108751b86443852b6741680d704168fTim Peters if a is None: 142d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Initialize from current time 143d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters import time 1440de88fc4b108751b86443852b6741680d704168fTim Peters a = long(time.time() * 256) 1450de88fc4b108751b86443852b6741680d704168fTim Peters 1460de88fc4b108751b86443852b6741680d704168fTim Peters if type(a) not in (type(3), type(3L)): 1470de88fc4b108751b86443852b6741680d704168fTim Peters a = hash(a) 1480de88fc4b108751b86443852b6741680d704168fTim Peters 1490de88fc4b108751b86443852b6741680d704168fTim Peters a, x = divmod(a, 30268) 1500de88fc4b108751b86443852b6741680d704168fTim Peters a, y = divmod(a, 30306) 1510de88fc4b108751b86443852b6741680d704168fTim Peters a, z = divmod(a, 30322) 1520de88fc4b108751b86443852b6741680d704168fTim Peters self._seed = int(x)+1, int(y)+1, int(z)+1 153d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 154cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def random(self): 155cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters """Get the next random number in the range [0.0, 1.0).""" 156cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 157cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Wichman-Hill random number generator. 158cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # 159cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Wichmann, B. A. & Hill, I. D. (1982) 160cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Algorithm AS 183: 161cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # An efficient and portable pseudo-random number generator 162cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Applied Statistics 31 (1982) 188-190 163cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # 164cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # see also: 165cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Correction to Algorithm AS 183 166cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Applied Statistics 33 (1984) 123 167cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # 168cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # McLeod, A. I. (1985) 169cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # A remark on Algorithm AS 183 170cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Applied Statistics 34 (1985),198-200 171cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 172cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # This part is thread-unsafe: 173cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # BEGIN CRITICAL SECTION 174cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters x, y, z = self._seed 175cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters x = (171 * x) % 30269 176cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters y = (172 * y) % 30307 177cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters z = (170 * z) % 30323 178cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters self._seed = x, y, z 179cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # END CRITICAL SECTION 180cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 181cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Note: on a platform using IEEE-754 double arithmetic, this can 182cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # never return 0.0 (asserted by Tim; proof too long for a comment). 183cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0 184cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 185d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def getstate(self): 186d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Return internal state; can be passed to setstate() later.""" 187d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return self.VERSION, self._seed, self.gauss_next 188d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 189d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def setstate(self, state): 190d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Restore internal state from object returned by getstate().""" 191d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters version = state[0] 192d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if version == 1: 193d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters version, self._seed, self.gauss_next = state 194d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 195d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError("state with version %s passed to " 196d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters "Random.setstate() of version %s" % 197d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters (version, self.VERSION)) 198d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 199d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters def jumpahead(self, n): 200d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters """Act as if n calls to random() were made, but quickly. 201d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters 202d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters n is an int, greater than or equal to 0. 203d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters 204d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters Example use: If you have 2 threads and know that each will 205d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters consume no more than a million random numbers, create two Random 206d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters objects r1 and r2, then do 207d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters r2.setstate(r1.getstate()) 208d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters r2.jumpahead(1000000) 209d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters Then r1 and r2 will use guaranteed-disjoint segments of the full 210d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters period. 211d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters """ 212d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters 213d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters if not n >= 0: 214d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters raise ValueError("n must be >= 0") 215d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters x, y, z = self._seed 216d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters x = int(x * pow(171, n, 30269)) % 30269 217d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters y = int(y * pow(172, n, 30307)) % 30307 218d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters z = int(z * pow(170, n, 30323)) % 30323 219d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters self._seed = x, y, z 220d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters 2210de88fc4b108751b86443852b6741680d704168fTim Peters def __whseed(self, x=0, y=0, z=0): 2220de88fc4b108751b86443852b6741680d704168fTim Peters """Set the Wichmann-Hill seed from (x, y, z). 2230de88fc4b108751b86443852b6741680d704168fTim Peters 2240de88fc4b108751b86443852b6741680d704168fTim Peters These must be integers in the range [0, 256). 2250de88fc4b108751b86443852b6741680d704168fTim Peters """ 2260de88fc4b108751b86443852b6741680d704168fTim Peters 2270de88fc4b108751b86443852b6741680d704168fTim Peters if not type(x) == type(y) == type(z) == type(0): 2280de88fc4b108751b86443852b6741680d704168fTim Peters raise TypeError('seeds must be integers') 2290de88fc4b108751b86443852b6741680d704168fTim Peters if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): 2300de88fc4b108751b86443852b6741680d704168fTim Peters raise ValueError('seeds must be in range(0, 256)') 2310de88fc4b108751b86443852b6741680d704168fTim Peters if 0 == x == y == z: 2320de88fc4b108751b86443852b6741680d704168fTim Peters # Initialize from current time 2330de88fc4b108751b86443852b6741680d704168fTim Peters import time 2340de88fc4b108751b86443852b6741680d704168fTim Peters t = long(time.time() * 256) 2350de88fc4b108751b86443852b6741680d704168fTim Peters t = int((t&0xffffff) ^ (t>>24)) 2360de88fc4b108751b86443852b6741680d704168fTim Peters t, x = divmod(t, 256) 2370de88fc4b108751b86443852b6741680d704168fTim Peters t, y = divmod(t, 256) 2380de88fc4b108751b86443852b6741680d704168fTim Peters t, z = divmod(t, 256) 2390de88fc4b108751b86443852b6741680d704168fTim Peters # Zero is a poor seed, so substitute 1 2400de88fc4b108751b86443852b6741680d704168fTim Peters self._seed = (x or 1, y or 1, z or 1) 2410de88fc4b108751b86443852b6741680d704168fTim Peters 2420de88fc4b108751b86443852b6741680d704168fTim Peters def whseed(self, a=None): 2430de88fc4b108751b86443852b6741680d704168fTim Peters """Seed from hashable object's hash code. 2440de88fc4b108751b86443852b6741680d704168fTim Peters 2450de88fc4b108751b86443852b6741680d704168fTim Peters None or no argument seeds from current time. It is not guaranteed 2460de88fc4b108751b86443852b6741680d704168fTim Peters that objects with distinct hash codes lead to distinct internal 2470de88fc4b108751b86443852b6741680d704168fTim Peters states. 2480de88fc4b108751b86443852b6741680d704168fTim Peters 2490de88fc4b108751b86443852b6741680d704168fTim Peters This is obsolete, provided for compatibility with the seed routine 2500de88fc4b108751b86443852b6741680d704168fTim Peters used prior to Python 2.1. Use the .seed() method instead. 2510de88fc4b108751b86443852b6741680d704168fTim Peters """ 2520de88fc4b108751b86443852b6741680d704168fTim Peters 2530de88fc4b108751b86443852b6741680d704168fTim Peters if a is None: 2540de88fc4b108751b86443852b6741680d704168fTim Peters self.__whseed() 2550de88fc4b108751b86443852b6741680d704168fTim Peters return 2560de88fc4b108751b86443852b6741680d704168fTim Peters a = hash(a) 2570de88fc4b108751b86443852b6741680d704168fTim Peters a, x = divmod(a, 256) 2580de88fc4b108751b86443852b6741680d704168fTim Peters a, y = divmod(a, 256) 2590de88fc4b108751b86443852b6741680d704168fTim Peters a, z = divmod(a, 256) 2600de88fc4b108751b86443852b6741680d704168fTim Peters x = (x + a) % 256 or 1 2610de88fc4b108751b86443852b6741680d704168fTim Peters y = (y + a) % 256 or 1 2620de88fc4b108751b86443852b6741680d704168fTim Peters z = (z + a) % 256 or 1 2630de88fc4b108751b86443852b6741680d704168fTim Peters self.__whseed(x, y, z) 2640de88fc4b108751b86443852b6741680d704168fTim Peters 265cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- Methods below this point do not need to be overridden when 266cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- subclassing for the purpose of using a different core generator. 267d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 268cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- pickle support ------------------- 269d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 270cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __getstate__(self): # for pickle 271cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters return self.getstate() 272d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 273cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __setstate__(self, state): # for pickle 274cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters self.setstate(state) 275cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 276cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- integer methods ------------------- 277d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 278d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def randrange(self, start, stop=None, step=1, int=int, default=None): 279d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random item from range(start, stop[, step]). 280d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 281d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters This fixes the problem with randint() which includes the 282d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters endpoint; in Python this is usually not what you want. 283d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Do not supply the 'int' and 'default' arguments. 284d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 285d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 286d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # This code is a bit messy to make it fast for the 287d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # common case while still doing adequate error checking 288d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istart = int(start) 289d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart != start: 290d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer arg 1 for randrange()" 291d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if stop is default: 292d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart > 0: 293d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return int(self.random() * istart) 294d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 295d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istop = int(stop) 296d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istop != stop: 297d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer stop for randrange()" 298d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if step == 1: 299d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart < istop: 300d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return istart + int(self.random() * 301d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters (istop - istart)) 302d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 303d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istep = int(step) 304d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep != step: 305d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer step for randrange()" 306d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep > 0: 307d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters n = (istop - istart + istep - 1) / istep 308d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif istep < 0: 309d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters n = (istop - istart + istep + 1) / istep 310d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 311d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "zero step for randrange()" 312d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 313d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if n <= 0: 314d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 315d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return istart + istep*int(self.random() * n) 316d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 317d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def randint(self, a, b): 318cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters """Return random integer in range [a, b], including both end points. 319d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 320cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters (Deprecated; use randrange(a, b+1).) 321d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 322d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 323d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return self.randrange(a, b+1) 324d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 325cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- sequence methods ------------------- 326cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 327d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def choice(self, seq): 328d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random element from a non-empty sequence.""" 329d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return seq[int(self.random() * len(seq))] 330d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 331d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def shuffle(self, x, random=None, int=int): 332d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """x, random=random.random -> shuffle list x in place; return None. 333d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 334d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional arg random is a 0-argument function returning a random 335d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters float in [0.0, 1.0); by default, the standard random.random. 336d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 337d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Note that for even rather small len(x), the total number of 338d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters permutations of x is larger than the period of most random number 339d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters generators; this implies that "most" permutations of a long 340d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters sequence can never be generated. 341d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 342d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 343d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if random is None: 344d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 345d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters for i in xrange(len(x)-1, 0, -1): 346cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # pick an element in x[:i+1] with which to exchange x[i] 347d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters j = int(random() * (i+1)) 348d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x[i], x[j] = x[j], x[i] 349d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 350cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- real-valued distributions ------------------- 351cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 352cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- uniform distribution ------------------- 353d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 354d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def uniform(self, a, b): 355d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Get a random number in the range [a, b).""" 356d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return a + (b-a) * self.random() 357ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 358cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- normal distribution -------------------- 359ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 360d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def normalvariate(self, mu, sigma): 361d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mu = mean, sigma = standard deviation 362d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 363d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses Kinderman and Monahan method. Reference: Kinderman, 364d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # A.J. and Monahan, J.F., "Computer generation of random 365d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables using the ratio of uniform deviates", ACM Trans 366d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Math Software, 3, (1977), pp257-260. 367d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 368d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 369d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while 1: 370d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 371d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u2 = random() 372d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = NV_MAGICCONST*(u1-0.5)/u2 373d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters zz = z*z/4.0 374d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if zz <= -_log(u2): 375d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 376d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 377ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 378cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- lognormal distribution -------------------- 379ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 380d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def lognormvariate(self, mu, sigma): 381d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return _exp(self.normalvariate(mu, sigma)) 382ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 383cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- circular uniform -------------------- 384ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 385d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def cunifvariate(self, mean, arc): 386d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mean: mean angle (in radians between 0 and pi) 387d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # arc: range of distribution (in radians between 0 and pi) 388ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 389d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return (mean + arc * (self.random() - 0.5)) % _pi 390ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 391cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- exponential distribution -------------------- 392ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 393d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def expovariate(self, lambd): 394d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lambd: rate lambd = 1/mean 395d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # ('lambda' is a Python reserved word) 396ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 397d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 3980c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 399d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 400d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 401d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return -_log(u)/lambd 402ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 403cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- von Mises distribution -------------------- 404ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 405d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def vonmisesvariate(self, mu, kappa): 406d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mu: mean angle (in radians between 0 and 2*pi) 407d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # kappa: concentration parameter kappa (>= 0) 408d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # if kappa = 0 generate uniform random angle 4095810297052003f28788f6790ac799fe8e5373494Guido van Rossum 410d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Based upon an algorithm published in: Fisher, N.I., 411d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # "Statistical Analysis of Circular Data", Cambridge 412d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # University Press, 1993. 4135810297052003f28788f6790ac799fe8e5373494Guido van Rossum 414d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Thanks to Magnus Kessler for a correction to the 415d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # implementation of step 4. 4165810297052003f28788f6790ac799fe8e5373494Guido van Rossum 417d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 418d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if kappa <= 1e-6: 419d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return TWOPI * random() 420ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 421d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa) 422d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (a - _sqrt(2.0 * a))/(2.0 * kappa) 423d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = (1.0 + b * b)/(2.0 * b) 424ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 425d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while 1: 426d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 427ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 428d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(_pi * u1) 429d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters f = (1.0 + r * z)/(r + z) 430d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters c = kappa * (r - f) 431ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 432d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u2 = random() 433ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 434d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)): 435d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 436ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 437d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u3 = random() 438d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if u3 > 0.5: 439d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters theta = (mu % TWOPI) + _acos(f) 440d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 441d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters theta = (mu % TWOPI) - _acos(f) 442ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 443d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return theta 444ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 445cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- gamma distribution -------------------- 446ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 447d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gammavariate(self, alpha, beta): 448d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # beta times standard gamma 449d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters ainv = _sqrt(2.0 * alpha - 1.0) 450d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return beta * self.stdgamma(alpha, ainv, alpha - LOG4, alpha + ainv) 451d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 452d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def stdgamma(self, alpha, ainv, bbb, ccc): 453d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # ainv = sqrt(2 * alpha - 1) 454d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # bbb = alpha - log(4) 455d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # ccc = alpha + ainv 456d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 457d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 458d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if alpha <= 0.0: 459d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, 'stdgamma: alpha must be > 0.0' 460d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 461d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if alpha > 1.0: 462d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 463d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses R.C.H. Cheng, "The generation of Gamma 464d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables with non-integral shape parameters", 465d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Applied Statistics, (1977), 26, No. 1, p71-74 466d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 467d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while 1: 468d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 469d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u2 = random() 470d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters v = _log(u1/(1.0-u1))/ainv 471d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = alpha*_exp(v) 472d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = u1*u1*u2 473d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = bbb+ccc*v-x 474d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 475d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return x 476d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 477d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif alpha == 1.0: 478d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # expovariate(1) 4790c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 480d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 481d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 482d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return -_log(u) 483d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 484d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: # alpha is between 0 and 1 (exclusive) 485d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 486d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 487d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 488d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while 1: 489d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 490d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (_e + alpha)/_e 491d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters p = b*u 492d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if p <= 1.0: 493d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = pow(p, 1.0/alpha) 494d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 495d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # p > 1 496d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = -_log((b-p)/alpha) 497d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 498d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not (((p <= 1.0) and (u1 > _exp(-x))) or 499d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters ((p > 1) and (u1 > pow(x, alpha - 1.0)))): 500d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 501d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return x 502ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 50395bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 504cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Gauss (faster alternative) -------------------- 50595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 506d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gauss(self, mu, sigma): 507d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 508d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # When x and y are two variables from [0, 1), uniformly 509d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # distributed, then 510d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 511d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # cos(2*pi*x)*sqrt(-2*log(1-y)) 512d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # sin(2*pi*x)*sqrt(-2*log(1-y)) 513d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 514d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # are two *independent* variables with normal distribution 515d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (mu = 0, sigma = 1). 516d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (Lambert Meertens) 517d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (corrected version; bug discovered by Mike Miller, fixed by LM) 518d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 519d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Multithreading note: When two threads call this function 520d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # simultaneously, it is possible that they will receive the 521d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # same return value. The window is very small though. To 522d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # avoid this, you have to use a lock around all calls. (I 523d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # didn't want to slow this down in the serial case by using a 524d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lock here.) 525d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 526d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 527d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = self.gauss_next 528d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = None 529d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if z is None: 530d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x2pi = random() * TWOPI 531d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters g2rad = _sqrt(-2.0 * _log(1.0 - random())) 532d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(x2pi) * g2rad 533d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = _sin(x2pi) * g2rad 534d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 535d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 53695bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 537cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- beta -------------------- 53885e2e4742d0a1accecd02058a7907df36308297eTim Peters## See 53985e2e4742d0a1accecd02058a7907df36308297eTim Peters## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470 54085e2e4742d0a1accecd02058a7907df36308297eTim Peters## for Ivan Frohne's insightful analysis of why the original implementation: 54185e2e4742d0a1accecd02058a7907df36308297eTim Peters## 54285e2e4742d0a1accecd02058a7907df36308297eTim Peters## def betavariate(self, alpha, beta): 54385e2e4742d0a1accecd02058a7907df36308297eTim Peters## # Discrete Event Simulation in C, pp 87-88. 54485e2e4742d0a1accecd02058a7907df36308297eTim Peters## 54585e2e4742d0a1accecd02058a7907df36308297eTim Peters## y = self.expovariate(alpha) 54685e2e4742d0a1accecd02058a7907df36308297eTim Peters## z = self.expovariate(1.0/beta) 54785e2e4742d0a1accecd02058a7907df36308297eTim Peters## return z/(y+z) 54885e2e4742d0a1accecd02058a7907df36308297eTim Peters## 54985e2e4742d0a1accecd02058a7907df36308297eTim Peters## was dead wrong, and how it probably got that way. 55095bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 551d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def betavariate(self, alpha, beta): 55285e2e4742d0a1accecd02058a7907df36308297eTim Peters # This version due to Janne Sinkkonen, and matches all the std 55385e2e4742d0a1accecd02058a7907df36308297eTim Peters # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 55485e2e4742d0a1accecd02058a7907df36308297eTim Peters y = self.gammavariate(alpha, 1.) 55585e2e4742d0a1accecd02058a7907df36308297eTim Peters if y == 0: 55685e2e4742d0a1accecd02058a7907df36308297eTim Peters return 0.0 55785e2e4742d0a1accecd02058a7907df36308297eTim Peters else: 55885e2e4742d0a1accecd02058a7907df36308297eTim Peters return y / (y + self.gammavariate(beta, 1.)) 55995bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 560cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Pareto -------------------- 561cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 562d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def paretovariate(self, alpha): 563d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 495 564cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 565d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = self.random() 566d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return 1.0 / pow(u, 1.0/alpha) 567cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 568cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Weibull -------------------- 569cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 570d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def weibullvariate(self, alpha, beta): 571d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 499; bug fix courtesy Bill Arms 572cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 573d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = self.random() 574d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return alpha * pow(-_log(u), 1.0/beta) 5756c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum 576cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- test program -------------------- 577ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 578d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersdef _test_generator(n, funccall): 5790c9886d589ddebf32de0ca3f027a173222ed383aTim Peters import time 5800c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print n, 'times', funccall 5810c9886d589ddebf32de0ca3f027a173222ed383aTim Peters code = compile(funccall, funccall, 'eval') 5820c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sum = 0.0 5830c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = 0.0 5840c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = 1e10 5850c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = -1e10 5860c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t0 = time.time() 5870c9886d589ddebf32de0ca3f027a173222ed383aTim Peters for i in range(n): 5880c9886d589ddebf32de0ca3f027a173222ed383aTim Peters x = eval(code) 5890c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sum = sum + x 5900c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = sqsum + x*x 5910c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = min(x, smallest) 5920c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = max(x, largest) 5930c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t1 = time.time() 5940c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print round(t1-t0, 3), 'sec,', 5950c9886d589ddebf32de0ca3f027a173222ed383aTim Peters avg = sum/n 596d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters stddev = _sqrt(sqsum/n - avg*avg) 5970c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print 'avg %g, stddev %g, min %g, max %g' % \ 5980c9886d589ddebf32de0ca3f027a173222ed383aTim Peters (avg, stddev, smallest, largest) 599ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 600d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersdef _test(N=200): 601d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters print 'TWOPI =', TWOPI 602d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters print 'LOG4 =', LOG4 603d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters print 'NV_MAGICCONST =', NV_MAGICCONST 604d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters print 'SG_MAGICCONST =', SG_MAGICCONST 605d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'random()') 606d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'normalvariate(0.0, 1.0)') 607d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'lognormvariate(0.0, 1.0)') 608d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'cunifvariate(0.0, 1.0)') 609d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'expovariate(1.0)') 610d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'vonmisesvariate(0.0, 1.0)') 611d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(0.5, 1.0)') 612d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(0.9, 1.0)') 613d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(1.0, 1.0)') 614d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(2.0, 1.0)') 615d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(20.0, 1.0)') 616d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(200.0, 1.0)') 617d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gauss(0.0, 1.0)') 618d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'betavariate(3.0, 3.0)') 619d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'paretovariate(1.0)') 620d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'weibullvariate(1.0, 1.0)') 621d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 622cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Test jumpahead. 623cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters s = getstate() 624cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters jumpahead(N) 625cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters r1 = random() 626cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # now do it the slow way 627cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters setstate(s) 628cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters for i in range(N): 629cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters random() 630cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters r2 = random() 631cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters if r1 != r2: 632cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters raise ValueError("jumpahead test failed " + `(N, r1, r2)`) 633cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 634715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# Create one instance, seeded from current time, and export its methods 635715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# as module-level functions. The functions are not threadsafe, and state 636715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# is shared across all uses (both in the user's code and in the Python 637715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# libraries), but that's fine for most programs and is easier for the 638715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# casual user than making them instantiate their own Random() instance. 639d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_inst = Random() 640d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersseed = _inst.seed 641d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandom = _inst.random 642d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersuniform = _inst.uniform 643d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandint = _inst.randint 644d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterschoice = _inst.choice 645d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandrange = _inst.randrange 646d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersshuffle = _inst.shuffle 647d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersnormalvariate = _inst.normalvariate 648d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterslognormvariate = _inst.lognormvariate 649d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterscunifvariate = _inst.cunifvariate 650d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersexpovariate = _inst.expovariate 651d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersvonmisesvariate = _inst.vonmisesvariate 652d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgammavariate = _inst.gammavariate 653d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersstdgamma = _inst.stdgamma 654d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgauss = _inst.gauss 655d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersbetavariate = _inst.betavariate 656d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersparetovariate = _inst.paretovariate 657d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersweibullvariate = _inst.weibullvariate 658d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgetstate = _inst.getstate 659d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterssetstate = _inst.setstate 660d52269bfd029c4a517ea74c17edd5c3a250c366cTim Petersjumpahead = _inst.jumpahead 6610de88fc4b108751b86443852b6741680d704168fTim Peterswhseed = _inst.whseed 662d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 663ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumif __name__ == '__main__': 664d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test() 665