random.py revision ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019
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"] 840e6d213177dd571dd634d286ae45c38eb06d63b9Tim Peters 85dc47a89ff19b932c1c794422a5223847b4e64f0eTim Petersdef _verify(name, computed, expected): 86d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if abs(computed - expected) > 1e-7: 87d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError( 88d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters "computed value for %s deviates too much " 89d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters "(computed %g, expected %g)" % (name, computed, expected)) 90ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 91d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersNV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) 92dc47a89ff19b932c1c794422a5223847b4e64f0eTim Peters_verify('NV_MAGICCONST', NV_MAGICCONST, 1.71552776992141) 9333d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 94d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersTWOPI = 2.0*_pi 95dc47a89ff19b932c1c794422a5223847b4e64f0eTim Peters_verify('TWOPI', TWOPI, 6.28318530718) 960c9886d589ddebf32de0ca3f027a173222ed383aTim Peters 97d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersLOG4 = _log(4.0) 98dc47a89ff19b932c1c794422a5223847b4e64f0eTim Peters_verify('LOG4', LOG4, 1.38629436111989) 990c9886d589ddebf32de0ca3f027a173222ed383aTim Peters 100d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersSG_MAGICCONST = 1.0 + _log(4.5) 101dc47a89ff19b932c1c794422a5223847b4e64f0eTim Peters_verify('SG_MAGICCONST', SG_MAGICCONST, 2.50407739677627) 10233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 103d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersdel _verify 10433d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 105d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters# Translated by Guido van Rossum from C source provided by 106d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters# Adrian Baddeley. 10733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 108d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersclass Random: 10933d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 110d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters VERSION = 1 # used by getstate/setstate 11133d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 112d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def __init__(self, x=None): 113d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Initialize an instance. 11433d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 115d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional argument x controls seeding, as for Random.seed(). 116d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 11733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 118d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.seed(x) 11933d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 120cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- core generator ------------------- 121cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 122d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Specific to Wichmann-Hill generator. Subclasses wishing to use a 123d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters # different core generator should override the seed(), random(), 124cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # getstate(), setstate() and jumpahead() methods. 125ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 1260de88fc4b108751b86443852b6741680d704168fTim Peters def seed(self, a=None): 1270de88fc4b108751b86443852b6741680d704168fTim Peters """Initialize internal state from hashable object. 128d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 1290de88fc4b108751b86443852b6741680d704168fTim Peters None or no argument seeds from current time. 1300de88fc4b108751b86443852b6741680d704168fTim Peters 131bcd725fc456faca13f4598f87c0517f917711cdaTim Peters If a is not None or an int or long, hash(a) is used instead. 1320de88fc4b108751b86443852b6741680d704168fTim Peters 1330de88fc4b108751b86443852b6741680d704168fTim Peters If a is an int or long, a is used directly. Distinct values between 1340de88fc4b108751b86443852b6741680d704168fTim Peters 0 and 27814431486575L inclusive are guaranteed to yield distinct 1350de88fc4b108751b86443852b6741680d704168fTim Peters internal states (this guarantee is specific to the default 1360de88fc4b108751b86443852b6741680d704168fTim Peters Wichmann-Hill generator). 137d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 138d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 1390de88fc4b108751b86443852b6741680d704168fTim Peters if a is None: 140d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Initialize from current time 141d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters import time 1420de88fc4b108751b86443852b6741680d704168fTim Peters a = long(time.time() * 256) 1430de88fc4b108751b86443852b6741680d704168fTim Peters 1440de88fc4b108751b86443852b6741680d704168fTim Peters if type(a) not in (type(3), type(3L)): 1450de88fc4b108751b86443852b6741680d704168fTim Peters a = hash(a) 1460de88fc4b108751b86443852b6741680d704168fTim Peters 1470de88fc4b108751b86443852b6741680d704168fTim Peters a, x = divmod(a, 30268) 1480de88fc4b108751b86443852b6741680d704168fTim Peters a, y = divmod(a, 30306) 1490de88fc4b108751b86443852b6741680d704168fTim Peters a, z = divmod(a, 30322) 1500de88fc4b108751b86443852b6741680d704168fTim Peters self._seed = int(x)+1, int(y)+1, int(z)+1 151d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 15246c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters self.gauss_next = None 15346c04e140cf26d1b44935c28c6f15ea467400d22Tim 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 24246c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters self.gauss_next = None 24346c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters 2440de88fc4b108751b86443852b6741680d704168fTim Peters def whseed(self, a=None): 2450de88fc4b108751b86443852b6741680d704168fTim Peters """Seed from hashable object's hash code. 2460de88fc4b108751b86443852b6741680d704168fTim Peters 2470de88fc4b108751b86443852b6741680d704168fTim Peters None or no argument seeds from current time. It is not guaranteed 2480de88fc4b108751b86443852b6741680d704168fTim Peters that objects with distinct hash codes lead to distinct internal 2490de88fc4b108751b86443852b6741680d704168fTim Peters states. 2500de88fc4b108751b86443852b6741680d704168fTim Peters 2510de88fc4b108751b86443852b6741680d704168fTim Peters This is obsolete, provided for compatibility with the seed routine 2520de88fc4b108751b86443852b6741680d704168fTim Peters used prior to Python 2.1. Use the .seed() method instead. 2530de88fc4b108751b86443852b6741680d704168fTim Peters """ 2540de88fc4b108751b86443852b6741680d704168fTim Peters 2550de88fc4b108751b86443852b6741680d704168fTim Peters if a is None: 2560de88fc4b108751b86443852b6741680d704168fTim Peters self.__whseed() 2570de88fc4b108751b86443852b6741680d704168fTim Peters return 2580de88fc4b108751b86443852b6741680d704168fTim Peters a = hash(a) 2590de88fc4b108751b86443852b6741680d704168fTim Peters a, x = divmod(a, 256) 2600de88fc4b108751b86443852b6741680d704168fTim Peters a, y = divmod(a, 256) 2610de88fc4b108751b86443852b6741680d704168fTim Peters a, z = divmod(a, 256) 2620de88fc4b108751b86443852b6741680d704168fTim Peters x = (x + a) % 256 or 1 2630de88fc4b108751b86443852b6741680d704168fTim Peters y = (y + a) % 256 or 1 2640de88fc4b108751b86443852b6741680d704168fTim Peters z = (z + a) % 256 or 1 2650de88fc4b108751b86443852b6741680d704168fTim Peters self.__whseed(x, y, z) 2660de88fc4b108751b86443852b6741680d704168fTim Peters 267cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- Methods below this point do not need to be overridden when 268cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- subclassing for the purpose of using a different core generator. 269d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 270cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- pickle support ------------------- 271d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 272cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __getstate__(self): # for pickle 273cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters return self.getstate() 274d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 275cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __setstate__(self, state): # for pickle 276cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters self.setstate(state) 277cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 278cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- integer methods ------------------- 279d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 280d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def randrange(self, start, stop=None, step=1, int=int, default=None): 281d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random item from range(start, stop[, step]). 282d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 283d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters This fixes the problem with randint() which includes the 284d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters endpoint; in Python this is usually not what you want. 285d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Do not supply the 'int' and 'default' arguments. 286d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 287d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 288d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # This code is a bit messy to make it fast for the 289d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # common case while still doing adequate error checking 290d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istart = int(start) 291d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart != start: 292d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer arg 1 for randrange()" 293d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if stop is default: 294d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart > 0: 295d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return int(self.random() * istart) 296d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 297d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istop = int(stop) 298d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istop != stop: 299d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer stop for randrange()" 300d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if step == 1: 301d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart < istop: 302d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return istart + int(self.random() * 303d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters (istop - istart)) 304d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 305d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istep = int(step) 306d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep != step: 307d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer step for randrange()" 308d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep > 0: 309d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters n = (istop - istart + istep - 1) / istep 310d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif istep < 0: 311d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters n = (istop - istart + istep + 1) / istep 312d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 313d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "zero step for randrange()" 314d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 315d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if n <= 0: 316d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 317d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return istart + istep*int(self.random() * n) 318d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 319d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def randint(self, a, b): 320cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters """Return random integer in range [a, b], including both end points. 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 449ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger return beta * self.stdgamma(alpha) 450d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 451ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger def stdgamma(self, alpha, *args): # *args for Py2.2 compatiblity 452d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 453d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if alpha <= 0.0: 454d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, 'stdgamma: alpha must be > 0.0' 455d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 456d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if alpha > 1.0: 457d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 458d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses R.C.H. Cheng, "The generation of Gamma 459d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables with non-integral shape parameters", 460d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Applied Statistics, (1977), 26, No. 1, p71-74 461d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 462ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger ainv = _sqrt(2.0 * alpha - 1.0) 463ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger bbb = alpha - LOG4 464ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger ccc = alpha + ainv 465ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger 466d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while 1: 467d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 468d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u2 = random() 469d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters v = _log(u1/(1.0-u1))/ainv 470d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = alpha*_exp(v) 471d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = u1*u1*u2 472d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = bbb+ccc*v-x 473d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 474d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return x 475d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 476d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif alpha == 1.0: 477d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # expovariate(1) 4780c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 479d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 480d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 481d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return -_log(u) 482d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 483d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: # alpha is between 0 and 1 (exclusive) 484d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 485d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 486d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 487d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while 1: 488d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 489d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (_e + alpha)/_e 490d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters p = b*u 491d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if p <= 1.0: 492d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = pow(p, 1.0/alpha) 493d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 494d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # p > 1 495d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = -_log((b-p)/alpha) 496d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 497d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not (((p <= 1.0) and (u1 > _exp(-x))) or 498d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters ((p > 1) and (u1 > pow(x, alpha - 1.0)))): 499d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 500d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return x 501ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 50295bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 503cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Gauss (faster alternative) -------------------- 50495bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 505d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gauss(self, mu, sigma): 506d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 507d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # When x and y are two variables from [0, 1), uniformly 508d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # distributed, then 509d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 510d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # cos(2*pi*x)*sqrt(-2*log(1-y)) 511d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # sin(2*pi*x)*sqrt(-2*log(1-y)) 512d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 513d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # are two *independent* variables with normal distribution 514d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (mu = 0, sigma = 1). 515d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (Lambert Meertens) 516d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (corrected version; bug discovered by Mike Miller, fixed by LM) 517d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 518d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Multithreading note: When two threads call this function 519d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # simultaneously, it is possible that they will receive the 520d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # same return value. The window is very small though. To 521d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # avoid this, you have to use a lock around all calls. (I 522d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # didn't want to slow this down in the serial case by using a 523d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lock here.) 524d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 525d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 526d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = self.gauss_next 527d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = None 528d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if z is None: 529d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x2pi = random() * TWOPI 530d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters g2rad = _sqrt(-2.0 * _log(1.0 - random())) 531d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(x2pi) * g2rad 532d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = _sin(x2pi) * g2rad 533d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 534d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 53595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 536cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- beta -------------------- 53785e2e4742d0a1accecd02058a7907df36308297eTim Peters## See 53885e2e4742d0a1accecd02058a7907df36308297eTim Peters## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470 53985e2e4742d0a1accecd02058a7907df36308297eTim Peters## for Ivan Frohne's insightful analysis of why the original implementation: 54085e2e4742d0a1accecd02058a7907df36308297eTim Peters## 54185e2e4742d0a1accecd02058a7907df36308297eTim Peters## def betavariate(self, alpha, beta): 54285e2e4742d0a1accecd02058a7907df36308297eTim Peters## # Discrete Event Simulation in C, pp 87-88. 54385e2e4742d0a1accecd02058a7907df36308297eTim Peters## 54485e2e4742d0a1accecd02058a7907df36308297eTim Peters## y = self.expovariate(alpha) 54585e2e4742d0a1accecd02058a7907df36308297eTim Peters## z = self.expovariate(1.0/beta) 54685e2e4742d0a1accecd02058a7907df36308297eTim Peters## return z/(y+z) 54785e2e4742d0a1accecd02058a7907df36308297eTim Peters## 54885e2e4742d0a1accecd02058a7907df36308297eTim Peters## was dead wrong, and how it probably got that way. 54995bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 550d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def betavariate(self, alpha, beta): 55185e2e4742d0a1accecd02058a7907df36308297eTim Peters # This version due to Janne Sinkkonen, and matches all the std 55285e2e4742d0a1accecd02058a7907df36308297eTim Peters # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 55385e2e4742d0a1accecd02058a7907df36308297eTim Peters y = self.gammavariate(alpha, 1.) 55485e2e4742d0a1accecd02058a7907df36308297eTim Peters if y == 0: 55585e2e4742d0a1accecd02058a7907df36308297eTim Peters return 0.0 55685e2e4742d0a1accecd02058a7907df36308297eTim Peters else: 55785e2e4742d0a1accecd02058a7907df36308297eTim Peters return y / (y + self.gammavariate(beta, 1.)) 55895bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 559cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Pareto -------------------- 560cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 561d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def paretovariate(self, alpha): 562d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 495 563cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 564d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = self.random() 565d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return 1.0 / pow(u, 1.0/alpha) 566cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 567cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Weibull -------------------- 568cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 569d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def weibullvariate(self, alpha, beta): 570d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 499; bug fix courtesy Bill Arms 571cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 572d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = self.random() 573d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return alpha * pow(-_log(u), 1.0/beta) 5746c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum 575cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- test program -------------------- 576ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 577d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersdef _test_generator(n, funccall): 5780c9886d589ddebf32de0ca3f027a173222ed383aTim Peters import time 5790c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print n, 'times', funccall 5800c9886d589ddebf32de0ca3f027a173222ed383aTim Peters code = compile(funccall, funccall, 'eval') 5810c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sum = 0.0 5820c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = 0.0 5830c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = 1e10 5840c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = -1e10 5850c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t0 = time.time() 5860c9886d589ddebf32de0ca3f027a173222ed383aTim Peters for i in range(n): 5870c9886d589ddebf32de0ca3f027a173222ed383aTim Peters x = eval(code) 5880c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sum = sum + x 5890c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = sqsum + x*x 5900c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = min(x, smallest) 5910c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = max(x, largest) 5920c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t1 = time.time() 5930c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print round(t1-t0, 3), 'sec,', 5940c9886d589ddebf32de0ca3f027a173222ed383aTim Peters avg = sum/n 595d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters stddev = _sqrt(sqsum/n - avg*avg) 5960c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print 'avg %g, stddev %g, min %g, max %g' % \ 5970c9886d589ddebf32de0ca3f027a173222ed383aTim Peters (avg, stddev, smallest, largest) 598ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 599d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersdef _test(N=200): 600d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters print 'TWOPI =', TWOPI 601d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters print 'LOG4 =', LOG4 602d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters print 'NV_MAGICCONST =', NV_MAGICCONST 603d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters print 'SG_MAGICCONST =', SG_MAGICCONST 604d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'random()') 605d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'normalvariate(0.0, 1.0)') 606d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'lognormvariate(0.0, 1.0)') 607d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'cunifvariate(0.0, 1.0)') 608d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'expovariate(1.0)') 609d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'vonmisesvariate(0.0, 1.0)') 610d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(0.5, 1.0)') 611d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(0.9, 1.0)') 612d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(1.0, 1.0)') 613d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(2.0, 1.0)') 614d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(20.0, 1.0)') 615d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(200.0, 1.0)') 616d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gauss(0.0, 1.0)') 617d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'betavariate(3.0, 3.0)') 618d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'paretovariate(1.0)') 619d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'weibullvariate(1.0, 1.0)') 620d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 621cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Test jumpahead. 622cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters s = getstate() 623cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters jumpahead(N) 624cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters r1 = random() 625cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # now do it the slow way 626cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters setstate(s) 627cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters for i in range(N): 628cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters random() 629cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters r2 = random() 630cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters if r1 != r2: 631cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters raise ValueError("jumpahead test failed " + `(N, r1, r2)`) 632cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 633715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# Create one instance, seeded from current time, and export its methods 634715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# as module-level functions. The functions are not threadsafe, and state 635715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# is shared across all uses (both in the user's code and in the Python 636715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# libraries), but that's fine for most programs and is easier for the 637715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# casual user than making them instantiate their own Random() instance. 638d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_inst = Random() 639d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersseed = _inst.seed 640d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandom = _inst.random 641d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersuniform = _inst.uniform 642d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandint = _inst.randint 643d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterschoice = _inst.choice 644d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandrange = _inst.randrange 645d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersshuffle = _inst.shuffle 646d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersnormalvariate = _inst.normalvariate 647d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterslognormvariate = _inst.lognormvariate 648d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterscunifvariate = _inst.cunifvariate 649d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersexpovariate = _inst.expovariate 650d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersvonmisesvariate = _inst.vonmisesvariate 651d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgammavariate = _inst.gammavariate 652d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersstdgamma = _inst.stdgamma 653d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgauss = _inst.gauss 654d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersbetavariate = _inst.betavariate 655d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersparetovariate = _inst.paretovariate 656d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersweibullvariate = _inst.weibullvariate 657d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgetstate = _inst.getstate 658d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterssetstate = _inst.setstate 659d52269bfd029c4a517ea74c17edd5c3a250c366cTim Petersjumpahead = _inst.jumpahead 6600de88fc4b108751b86443852b6741680d704168fTim Peterswhseed = _inst.whseed 661d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 662ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumif __name__ == '__main__': 663d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test() 664