random.py revision ef4d4bdc3c99d3120a401b7af6e06610716d2e47
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: 109c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Random number generator base class used by bound module functions. 110c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 111c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Used to instantiate instances of Random to get generators that don't 112c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger share state. Especially useful for multi-threaded programs, creating 113c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger a different instance of Random for each thread, and using the jumpahead() 114c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger method to ensure that the generated sequences seen by each thread don't 115c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger overlap. 116c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 117c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Class Random can also be subclassed if you want to use a different basic 118c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger generator of your own devising: in that case, override the following 119c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger methods: random(), seed(), getstate(), setstate() and jumpahead(). 120ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 121c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 12233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 123d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters VERSION = 1 # used by getstate/setstate 12433d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 125d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def __init__(self, x=None): 126d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Initialize an instance. 12733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 128d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional argument x controls seeding, as for Random.seed(). 129d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 13033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 131d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.seed(x) 13233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 133cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- core generator ------------------- 134cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 135d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Specific to Wichmann-Hill generator. Subclasses wishing to use a 136d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters # different core generator should override the seed(), random(), 137cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # getstate(), setstate() and jumpahead() methods. 138ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 1390de88fc4b108751b86443852b6741680d704168fTim Peters def seed(self, a=None): 1400de88fc4b108751b86443852b6741680d704168fTim Peters """Initialize internal state from hashable object. 141d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 1420de88fc4b108751b86443852b6741680d704168fTim Peters None or no argument seeds from current time. 1430de88fc4b108751b86443852b6741680d704168fTim Peters 144bcd725fc456faca13f4598f87c0517f917711cdaTim Peters If a is not None or an int or long, hash(a) is used instead. 1450de88fc4b108751b86443852b6741680d704168fTim Peters 1460de88fc4b108751b86443852b6741680d704168fTim Peters If a is an int or long, a is used directly. Distinct values between 1470de88fc4b108751b86443852b6741680d704168fTim Peters 0 and 27814431486575L inclusive are guaranteed to yield distinct 1480de88fc4b108751b86443852b6741680d704168fTim Peters internal states (this guarantee is specific to the default 1490de88fc4b108751b86443852b6741680d704168fTim Peters Wichmann-Hill generator). 150d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 151d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 1520de88fc4b108751b86443852b6741680d704168fTim Peters if a is None: 153d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Initialize from current time 154d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters import time 1550de88fc4b108751b86443852b6741680d704168fTim Peters a = long(time.time() * 256) 1560de88fc4b108751b86443852b6741680d704168fTim Peters 1570de88fc4b108751b86443852b6741680d704168fTim Peters if type(a) not in (type(3), type(3L)): 1580de88fc4b108751b86443852b6741680d704168fTim Peters a = hash(a) 1590de88fc4b108751b86443852b6741680d704168fTim Peters 1600de88fc4b108751b86443852b6741680d704168fTim Peters a, x = divmod(a, 30268) 1610de88fc4b108751b86443852b6741680d704168fTim Peters a, y = divmod(a, 30306) 1620de88fc4b108751b86443852b6741680d704168fTim Peters a, z = divmod(a, 30322) 1630de88fc4b108751b86443852b6741680d704168fTim Peters self._seed = int(x)+1, int(y)+1, int(z)+1 164d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 16546c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters self.gauss_next = None 16646c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters 167cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def random(self): 168cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters """Get the next random number in the range [0.0, 1.0).""" 169cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 170cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Wichman-Hill random number generator. 171cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # 172cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Wichmann, B. A. & Hill, I. D. (1982) 173cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Algorithm AS 183: 174cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # An efficient and portable pseudo-random number generator 175cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Applied Statistics 31 (1982) 188-190 176cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # 177cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # see also: 178cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Correction to Algorithm AS 183 179cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Applied Statistics 33 (1984) 123 180cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # 181cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # McLeod, A. I. (1985) 182cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # A remark on Algorithm AS 183 183cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Applied Statistics 34 (1985),198-200 184cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 185cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # This part is thread-unsafe: 186cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # BEGIN CRITICAL SECTION 187cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters x, y, z = self._seed 188cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters x = (171 * x) % 30269 189cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters y = (172 * y) % 30307 190cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters z = (170 * z) % 30323 191cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters self._seed = x, y, z 192cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # END CRITICAL SECTION 193cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 194cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Note: on a platform using IEEE-754 double arithmetic, this can 195cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # never return 0.0 (asserted by Tim; proof too long for a comment). 196cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0 197cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 198d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def getstate(self): 199d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Return internal state; can be passed to setstate() later.""" 200d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return self.VERSION, self._seed, self.gauss_next 201d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 202d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def setstate(self, state): 203d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Restore internal state from object returned by getstate().""" 204d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters version = state[0] 205d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if version == 1: 206d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters version, self._seed, self.gauss_next = state 207d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 208d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError("state with version %s passed to " 209d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters "Random.setstate() of version %s" % 210d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters (version, self.VERSION)) 211d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 212d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters def jumpahead(self, n): 213d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters """Act as if n calls to random() were made, but quickly. 214d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters 215d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters n is an int, greater than or equal to 0. 216d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters 217d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters Example use: If you have 2 threads and know that each will 218d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters consume no more than a million random numbers, create two Random 219d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters objects r1 and r2, then do 220d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters r2.setstate(r1.getstate()) 221d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters r2.jumpahead(1000000) 222d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters Then r1 and r2 will use guaranteed-disjoint segments of the full 223d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters period. 224d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters """ 225d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters 226d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters if not n >= 0: 227d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters raise ValueError("n must be >= 0") 228d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters x, y, z = self._seed 229d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters x = int(x * pow(171, n, 30269)) % 30269 230d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters y = int(y * pow(172, n, 30307)) % 30307 231d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters z = int(z * pow(170, n, 30323)) % 30323 232d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters self._seed = x, y, z 233d52269bfd029c4a517ea74c17edd5c3a250c366cTim Peters 2340de88fc4b108751b86443852b6741680d704168fTim Peters def __whseed(self, x=0, y=0, z=0): 2350de88fc4b108751b86443852b6741680d704168fTim Peters """Set the Wichmann-Hill seed from (x, y, z). 2360de88fc4b108751b86443852b6741680d704168fTim Peters 2370de88fc4b108751b86443852b6741680d704168fTim Peters These must be integers in the range [0, 256). 2380de88fc4b108751b86443852b6741680d704168fTim Peters """ 2390de88fc4b108751b86443852b6741680d704168fTim Peters 2400de88fc4b108751b86443852b6741680d704168fTim Peters if not type(x) == type(y) == type(z) == type(0): 2410de88fc4b108751b86443852b6741680d704168fTim Peters raise TypeError('seeds must be integers') 2420de88fc4b108751b86443852b6741680d704168fTim Peters if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): 2430de88fc4b108751b86443852b6741680d704168fTim Peters raise ValueError('seeds must be in range(0, 256)') 2440de88fc4b108751b86443852b6741680d704168fTim Peters if 0 == x == y == z: 2450de88fc4b108751b86443852b6741680d704168fTim Peters # Initialize from current time 2460de88fc4b108751b86443852b6741680d704168fTim Peters import time 2470de88fc4b108751b86443852b6741680d704168fTim Peters t = long(time.time() * 256) 2480de88fc4b108751b86443852b6741680d704168fTim Peters t = int((t&0xffffff) ^ (t>>24)) 2490de88fc4b108751b86443852b6741680d704168fTim Peters t, x = divmod(t, 256) 2500de88fc4b108751b86443852b6741680d704168fTim Peters t, y = divmod(t, 256) 2510de88fc4b108751b86443852b6741680d704168fTim Peters t, z = divmod(t, 256) 2520de88fc4b108751b86443852b6741680d704168fTim Peters # Zero is a poor seed, so substitute 1 2530de88fc4b108751b86443852b6741680d704168fTim Peters self._seed = (x or 1, y or 1, z or 1) 2540de88fc4b108751b86443852b6741680d704168fTim Peters 25546c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters self.gauss_next = None 25646c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters 2570de88fc4b108751b86443852b6741680d704168fTim Peters def whseed(self, a=None): 2580de88fc4b108751b86443852b6741680d704168fTim Peters """Seed from hashable object's hash code. 2590de88fc4b108751b86443852b6741680d704168fTim Peters 2600de88fc4b108751b86443852b6741680d704168fTim Peters None or no argument seeds from current time. It is not guaranteed 2610de88fc4b108751b86443852b6741680d704168fTim Peters that objects with distinct hash codes lead to distinct internal 2620de88fc4b108751b86443852b6741680d704168fTim Peters states. 2630de88fc4b108751b86443852b6741680d704168fTim Peters 2640de88fc4b108751b86443852b6741680d704168fTim Peters This is obsolete, provided for compatibility with the seed routine 2650de88fc4b108751b86443852b6741680d704168fTim Peters used prior to Python 2.1. Use the .seed() method instead. 2660de88fc4b108751b86443852b6741680d704168fTim Peters """ 2670de88fc4b108751b86443852b6741680d704168fTim Peters 2680de88fc4b108751b86443852b6741680d704168fTim Peters if a is None: 2690de88fc4b108751b86443852b6741680d704168fTim Peters self.__whseed() 2700de88fc4b108751b86443852b6741680d704168fTim Peters return 2710de88fc4b108751b86443852b6741680d704168fTim Peters a = hash(a) 2720de88fc4b108751b86443852b6741680d704168fTim Peters a, x = divmod(a, 256) 2730de88fc4b108751b86443852b6741680d704168fTim Peters a, y = divmod(a, 256) 2740de88fc4b108751b86443852b6741680d704168fTim Peters a, z = divmod(a, 256) 2750de88fc4b108751b86443852b6741680d704168fTim Peters x = (x + a) % 256 or 1 2760de88fc4b108751b86443852b6741680d704168fTim Peters y = (y + a) % 256 or 1 2770de88fc4b108751b86443852b6741680d704168fTim Peters z = (z + a) % 256 or 1 2780de88fc4b108751b86443852b6741680d704168fTim Peters self.__whseed(x, y, z) 2790de88fc4b108751b86443852b6741680d704168fTim Peters 280cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- Methods below this point do not need to be overridden when 281cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- subclassing for the purpose of using a different core generator. 282d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 283cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- pickle support ------------------- 284d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 285cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __getstate__(self): # for pickle 286cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters return self.getstate() 287d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 288cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __setstate__(self, state): # for pickle 289cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters self.setstate(state) 290cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 291cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- integer methods ------------------- 292d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 293d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def randrange(self, start, stop=None, step=1, int=int, default=None): 294d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random item from range(start, stop[, step]). 295d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 296d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters This fixes the problem with randint() which includes the 297d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters endpoint; in Python this is usually not what you want. 298d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Do not supply the 'int' and 'default' arguments. 299d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 300d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 301d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # This code is a bit messy to make it fast for the 302d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # common case while still doing adequate error checking 303d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istart = int(start) 304d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart != start: 305d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer arg 1 for randrange()" 306d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if stop is default: 307d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart > 0: 308d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return int(self.random() * istart) 309d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 310d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istop = int(stop) 311d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istop != stop: 312d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer stop for randrange()" 313d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if step == 1: 314d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart < istop: 315d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return istart + int(self.random() * 316d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters (istop - istart)) 317d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 318d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istep = int(step) 319d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep != step: 320d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer step for randrange()" 321d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep > 0: 322d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters n = (istop - istart + istep - 1) / istep 323d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif istep < 0: 324d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters n = (istop - istart + istep + 1) / istep 325d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 326d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "zero step for randrange()" 327d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 328d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if n <= 0: 329d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 330d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return istart + istep*int(self.random() * n) 331d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 332d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def randint(self, a, b): 333cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters """Return random integer in range [a, b], including both end points. 334d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 335d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 336d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return self.randrange(a, b+1) 337d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 338cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- sequence methods ------------------- 339cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 340d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def choice(self, seq): 341d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random element from a non-empty sequence.""" 342d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return seq[int(self.random() * len(seq))] 343d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 344d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def shuffle(self, x, random=None, int=int): 345d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """x, random=random.random -> shuffle list x in place; return None. 346d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 347d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional arg random is a 0-argument function returning a random 348d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters float in [0.0, 1.0); by default, the standard random.random. 349d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 350d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Note that for even rather small len(x), the total number of 351d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters permutations of x is larger than the period of most random number 352d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters generators; this implies that "most" permutations of a long 353d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters sequence can never be generated. 354d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 355d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 356d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if random is None: 357d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 358d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters for i in xrange(len(x)-1, 0, -1): 359cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # pick an element in x[:i+1] with which to exchange x[i] 360d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters j = int(random() * (i+1)) 361d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x[i], x[j] = x[j], x[i] 362d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 363cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- real-valued distributions ------------------- 364cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 365cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- uniform distribution ------------------- 366d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 367d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def uniform(self, a, b): 368d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Get a random number in the range [a, b).""" 369d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return a + (b-a) * self.random() 370ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 371cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- normal distribution -------------------- 372ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 373d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def normalvariate(self, mu, sigma): 374c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Normal distribution. 375c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 376c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean, and sigma is the standard deviation. 377ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 378c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 379d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mu = mean, sigma = standard deviation 380d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 381d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses Kinderman and Monahan method. Reference: Kinderman, 382d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # A.J. and Monahan, J.F., "Computer generation of random 383d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables using the ratio of uniform deviates", ACM Trans 384d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Math Software, 3, (1977), pp257-260. 385d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 386d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 387d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while 1: 388d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 389d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u2 = random() 390d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = NV_MAGICCONST*(u1-0.5)/u2 391d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters zz = z*z/4.0 392d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if zz <= -_log(u2): 393d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 394d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 395ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 396cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- lognormal distribution -------------------- 397ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 398d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def lognormvariate(self, mu, sigma): 399c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Log normal distribution. 400c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 401c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger If you take the natural logarithm of this distribution, you'll get a 402c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger normal distribution with mean mu and standard deviation sigma. 403c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu can have any value, and sigma must be greater than zero. 404ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 405c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 406d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return _exp(self.normalvariate(mu, sigma)) 407ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 408cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- circular uniform -------------------- 409ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 410d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def cunifvariate(self, mean, arc): 411c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Circular uniform distribution. 412c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 413c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mean is the mean angle, and arc is the range of the distribution, 414c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger centered around the mean angle. Both values must be expressed in 415c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger radians. Returned values range between mean - arc/2 and 416c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mean + arc/2 and are normalized to between 0 and pi. 417c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 418c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Deprecated in version 2.3. Use: 419c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger (mean + arc * (Random.random() - 0.5)) % Math.pi 420ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 421c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 422d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mean: mean angle (in radians between 0 and pi) 423d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # arc: range of distribution (in radians between 0 and pi) 424c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger import warnings 425c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger warnings.warn("The cunifvariate function is deprecated; Use (mean " 426c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger "+ arc * (Random.random() - 0.5)) % Math.pi instead", 427c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger DeprecationWarning) 428ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 429d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return (mean + arc * (self.random() - 0.5)) % _pi 430ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 431cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- exponential distribution -------------------- 432ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 433d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def expovariate(self, lambd): 434c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Exponential distribution. 435c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 436c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger lambd is 1.0 divided by the desired mean. (The parameter would be 437c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger called "lambda", but that is a reserved word in Python.) Returned 438c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger values range from 0 to positive infinity. 439ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 440c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 441d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lambd: rate lambd = 1/mean 442d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # ('lambda' is a Python reserved word) 443ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 444d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 4450c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 446d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 447d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 448d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return -_log(u)/lambd 449ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 450cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- von Mises distribution -------------------- 451ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 452d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def vonmisesvariate(self, mu, kappa): 453c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Circular data distribution. 454ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 455c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean angle, expressed in radians between 0 and 2*pi, and 456c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger kappa is the concentration parameter, which must be greater than or 457c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger equal to zero. If kappa is equal to zero, this distribution reduces 458c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger to a uniform random angle over the range 0 to 2*pi. 459ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 460c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 461d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mu: mean angle (in radians between 0 and 2*pi) 462d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # kappa: concentration parameter kappa (>= 0) 463d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # if kappa = 0 generate uniform random angle 4645810297052003f28788f6790ac799fe8e5373494Guido van Rossum 465d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Based upon an algorithm published in: Fisher, N.I., 466d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # "Statistical Analysis of Circular Data", Cambridge 467d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # University Press, 1993. 4685810297052003f28788f6790ac799fe8e5373494Guido van Rossum 469d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Thanks to Magnus Kessler for a correction to the 470d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # implementation of step 4. 4715810297052003f28788f6790ac799fe8e5373494Guido van Rossum 472d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 473d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if kappa <= 1e-6: 474d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return TWOPI * random() 475ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 476d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa) 477d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (a - _sqrt(2.0 * a))/(2.0 * kappa) 478d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = (1.0 + b * b)/(2.0 * b) 479ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 480d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while 1: 481d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 482ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 483d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(_pi * u1) 484d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters f = (1.0 + r * z)/(r + z) 485d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters c = kappa * (r - f) 486ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 487d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u2 = random() 488ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 489d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)): 490d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 491ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 492d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u3 = random() 493d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if u3 > 0.5: 494d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters theta = (mu % TWOPI) + _acos(f) 495d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 496d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters theta = (mu % TWOPI) - _acos(f) 497ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 498d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return theta 499ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 500cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- gamma distribution -------------------- 501ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 502d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gammavariate(self, alpha, beta): 503c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Gamma distribution. Not the gamma function! 504c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 505c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Conditions on the parameters are alpha > 0 and beta > 0. 506c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 507c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 5088ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 509b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 5108ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 511570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum # Warning: a few older sources define the gamma distribution in terms 512570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum # of alpha > -1.0 513570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum if alpha <= 0.0 or beta <= 0.0: 514570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum raise ValueError, 'gammavariate: alpha and beta must be > 0.0' 5158ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 516d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 517d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if alpha > 1.0: 518d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 519d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses R.C.H. Cheng, "The generation of Gamma 520d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables with non-integral shape parameters", 521d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Applied Statistics, (1977), 26, No. 1, p71-74 522d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 523ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger ainv = _sqrt(2.0 * alpha - 1.0) 524ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger bbb = alpha - LOG4 525ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger ccc = alpha + ainv 5268ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 527d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while 1: 528d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 529d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u2 = random() 530d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters v = _log(u1/(1.0-u1))/ainv 531d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = alpha*_exp(v) 532d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = u1*u1*u2 533d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = bbb+ccc*v-x 534d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 535b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return x * beta 536d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 537d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif alpha == 1.0: 538d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # expovariate(1) 5390c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 540d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 541d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 542b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return -_log(u) * beta 543d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 544d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: # alpha is between 0 and 1 (exclusive) 545d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 546d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 547d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 548d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while 1: 549d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 550d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (_e + alpha)/_e 551d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters p = b*u 552d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if p <= 1.0: 553d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = pow(p, 1.0/alpha) 554d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 555d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # p > 1 556d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = -_log((b-p)/alpha) 557d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 558d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not (((p <= 1.0) and (u1 > _exp(-x))) or 559d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters ((p > 1) and (u1 > pow(x, alpha - 1.0)))): 560d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 561b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return x * beta 562b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger 563b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger 564b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger def stdgamma(self, alpha, ainv, bbb, ccc): 565b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # This method was (and shall remain) undocumented. 566b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # This method is deprecated 567b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # for the following reasons: 568b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # 1. Returns same as .gammavariate(alpha, 1.0) 569b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # 2. Requires caller to provide 3 extra arguments 570b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # that are functions of alpha anyway 571b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # 3. Can't be used for alpha < 0.5 572b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger 573b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # ainv = sqrt(2 * alpha - 1) 574b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # bbb = alpha - log(4) 575b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # ccc = alpha + ainv 576b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger import warnings 577b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger warnings.warn("The stdgamma function is deprecated; " 578b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger "use gammavariate() instead", 579b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger DeprecationWarning) 580b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return self.gammavariate(alpha, 1.0) 581b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger 582ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 58395bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 584cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Gauss (faster alternative) -------------------- 58595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 586d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gauss(self, mu, sigma): 587c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Gaussian distribution. 588c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 589c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean, and sigma is the standard deviation. This is 590c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger slightly faster than the normalvariate() function. 591c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 592c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Not thread-safe without a lock around calls. 593ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 594c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 595d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 596d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # When x and y are two variables from [0, 1), uniformly 597d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # distributed, then 598d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 599d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # cos(2*pi*x)*sqrt(-2*log(1-y)) 600d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # sin(2*pi*x)*sqrt(-2*log(1-y)) 601d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 602d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # are two *independent* variables with normal distribution 603d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (mu = 0, sigma = 1). 604d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (Lambert Meertens) 605d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (corrected version; bug discovered by Mike Miller, fixed by LM) 606d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 607d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Multithreading note: When two threads call this function 608d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # simultaneously, it is possible that they will receive the 609d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # same return value. The window is very small though. To 610d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # avoid this, you have to use a lock around all calls. (I 611d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # didn't want to slow this down in the serial case by using a 612d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lock here.) 613d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 614d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 615d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = self.gauss_next 616d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = None 617d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if z is None: 618d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x2pi = random() * TWOPI 619d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters g2rad = _sqrt(-2.0 * _log(1.0 - random())) 620d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(x2pi) * g2rad 621d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = _sin(x2pi) * g2rad 622d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 623d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 62495bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 625cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- beta -------------------- 62685e2e4742d0a1accecd02058a7907df36308297eTim Peters## See 62785e2e4742d0a1accecd02058a7907df36308297eTim Peters## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470 62885e2e4742d0a1accecd02058a7907df36308297eTim Peters## for Ivan Frohne's insightful analysis of why the original implementation: 62985e2e4742d0a1accecd02058a7907df36308297eTim Peters## 63085e2e4742d0a1accecd02058a7907df36308297eTim Peters## def betavariate(self, alpha, beta): 63185e2e4742d0a1accecd02058a7907df36308297eTim Peters## # Discrete Event Simulation in C, pp 87-88. 63285e2e4742d0a1accecd02058a7907df36308297eTim Peters## 63385e2e4742d0a1accecd02058a7907df36308297eTim Peters## y = self.expovariate(alpha) 63485e2e4742d0a1accecd02058a7907df36308297eTim Peters## z = self.expovariate(1.0/beta) 63585e2e4742d0a1accecd02058a7907df36308297eTim Peters## return z/(y+z) 63685e2e4742d0a1accecd02058a7907df36308297eTim Peters## 63785e2e4742d0a1accecd02058a7907df36308297eTim Peters## was dead wrong, and how it probably got that way. 63895bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 639d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def betavariate(self, alpha, beta): 640c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Beta distribution. 641c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 642c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Conditions on the parameters are alpha > -1 and beta} > -1. 643c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Returned values range between 0 and 1. 644ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 645c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 646ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 64785e2e4742d0a1accecd02058a7907df36308297eTim Peters # This version due to Janne Sinkkonen, and matches all the std 64885e2e4742d0a1accecd02058a7907df36308297eTim Peters # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 64985e2e4742d0a1accecd02058a7907df36308297eTim Peters y = self.gammavariate(alpha, 1.) 65085e2e4742d0a1accecd02058a7907df36308297eTim Peters if y == 0: 65185e2e4742d0a1accecd02058a7907df36308297eTim Peters return 0.0 65285e2e4742d0a1accecd02058a7907df36308297eTim Peters else: 65385e2e4742d0a1accecd02058a7907df36308297eTim Peters return y / (y + self.gammavariate(beta, 1.)) 65495bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 655cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Pareto -------------------- 656cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 657d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def paretovariate(self, alpha): 658c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Pareto distribution. alpha is the shape parameter.""" 659d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 495 660cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 661d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = self.random() 662d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return 1.0 / pow(u, 1.0/alpha) 663cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 664cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Weibull -------------------- 665cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 666d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def weibullvariate(self, alpha, beta): 667c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Weibull distribution. 668c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 669c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger alpha is the scale parameter and beta is the shape parameter. 670ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 671c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 672d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 499; bug fix courtesy Bill Arms 673cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 674d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = self.random() 675d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return alpha * pow(-_log(u), 1.0/beta) 6766c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum 677cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- test program -------------------- 678ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 679d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersdef _test_generator(n, funccall): 6800c9886d589ddebf32de0ca3f027a173222ed383aTim Peters import time 6810c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print n, 'times', funccall 6820c9886d589ddebf32de0ca3f027a173222ed383aTim Peters code = compile(funccall, funccall, 'eval') 6830c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sum = 0.0 6840c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = 0.0 6850c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = 1e10 6860c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = -1e10 6870c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t0 = time.time() 6880c9886d589ddebf32de0ca3f027a173222ed383aTim Peters for i in range(n): 6890c9886d589ddebf32de0ca3f027a173222ed383aTim Peters x = eval(code) 6900c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sum = sum + x 6910c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = sqsum + x*x 6920c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = min(x, smallest) 6930c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = max(x, largest) 6940c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t1 = time.time() 6950c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print round(t1-t0, 3), 'sec,', 6960c9886d589ddebf32de0ca3f027a173222ed383aTim Peters avg = sum/n 697d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters stddev = _sqrt(sqsum/n - avg*avg) 6980c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print 'avg %g, stddev %g, min %g, max %g' % \ 6990c9886d589ddebf32de0ca3f027a173222ed383aTim Peters (avg, stddev, smallest, largest) 700ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 701b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettingerdef _test(N=20000): 702d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters print 'TWOPI =', TWOPI 703d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters print 'LOG4 =', LOG4 704d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters print 'NV_MAGICCONST =', NV_MAGICCONST 705d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters print 'SG_MAGICCONST =', SG_MAGICCONST 706d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'random()') 707d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'normalvariate(0.0, 1.0)') 708d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'lognormvariate(0.0, 1.0)') 709d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'cunifvariate(0.0, 1.0)') 710d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'expovariate(1.0)') 711d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'vonmisesvariate(0.0, 1.0)') 712b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger _test_generator(N, 'gammavariate(0.01, 1.0)') 713b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger _test_generator(N, 'gammavariate(0.1, 1.0)') 7148ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters _test_generator(N, 'gammavariate(0.1, 2.0)') 715d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(0.5, 1.0)') 716d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(0.9, 1.0)') 717d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(1.0, 1.0)') 718d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(2.0, 1.0)') 719d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(20.0, 1.0)') 720d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(200.0, 1.0)') 721d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gauss(0.0, 1.0)') 722d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'betavariate(3.0, 3.0)') 723d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'paretovariate(1.0)') 724d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'weibullvariate(1.0, 1.0)') 725d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 726cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # Test jumpahead. 727cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters s = getstate() 728cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters jumpahead(N) 729cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters r1 = random() 730cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # now do it the slow way 731cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters setstate(s) 732cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters for i in range(N): 733cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters random() 734cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters r2 = random() 735cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters if r1 != r2: 736cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters raise ValueError("jumpahead test failed " + `(N, r1, r2)`) 737cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 738715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# Create one instance, seeded from current time, and export its methods 739715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# as module-level functions. The functions are not threadsafe, and state 740715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# is shared across all uses (both in the user's code and in the Python 741715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# libraries), but that's fine for most programs and is easier for the 742715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# casual user than making them instantiate their own Random() instance. 743d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_inst = Random() 744d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersseed = _inst.seed 745d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandom = _inst.random 746d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersuniform = _inst.uniform 747d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandint = _inst.randint 748d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterschoice = _inst.choice 749d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandrange = _inst.randrange 750d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersshuffle = _inst.shuffle 751d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersnormalvariate = _inst.normalvariate 752d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterslognormvariate = _inst.lognormvariate 753d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterscunifvariate = _inst.cunifvariate 754d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersexpovariate = _inst.expovariate 755d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersvonmisesvariate = _inst.vonmisesvariate 756d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgammavariate = _inst.gammavariate 757d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersstdgamma = _inst.stdgamma 758d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgauss = _inst.gauss 759d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersbetavariate = _inst.betavariate 760d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersparetovariate = _inst.paretovariate 761d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersweibullvariate = _inst.weibullvariate 762d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgetstate = _inst.getstate 763d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterssetstate = _inst.setstate 764d52269bfd029c4a517ea74c17edd5c3a250c366cTim Petersjumpahead = _inst.jumpahead 7650de88fc4b108751b86443852b6741680d704168fTim Peterswhseed = _inst.whseed 766d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 767ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumif __name__ == '__main__': 768d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test() 769