+++ /dev/null
-"""Random variable generators.\r
-\r
- integers\r
- --------\r
- uniform within range\r
-\r
- sequences\r
- ---------\r
- pick random element\r
- pick random sample\r
- generate random permutation\r
-\r
- distributions on the real line:\r
- ------------------------------\r
- uniform\r
- triangular\r
- normal (Gaussian)\r
- lognormal\r
- negative exponential\r
- gamma\r
- beta\r
- pareto\r
- Weibull\r
-\r
- distributions on the circle (angles 0 to 2pi)\r
- ---------------------------------------------\r
- circular uniform\r
- von Mises\r
-\r
-General notes on the underlying Mersenne Twister core generator:\r
-\r
-* The period is 2**19937-1.\r
-* It is one of the most extensively tested generators in existence.\r
-* Without a direct way to compute N steps forward, the semantics of\r
- jumpahead(n) are weakened to simply jump to another distant state and rely\r
- on the large period to avoid overlapping sequences.\r
-* The random() method is implemented in C, executes in a single Python step,\r
- and is, therefore, threadsafe.\r
-\r
-"""\r
-\r
-from __future__ import division\r
-from warnings import warn as _warn\r
-from types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType\r
-from math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil\r
-from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin\r
-from os import urandom as _urandom\r
-from binascii import hexlify as _hexlify\r
-import hashlib as _hashlib\r
-\r
-__all__ = ["Random","seed","random","uniform","randint","choice","sample",\r
- "randrange","shuffle","normalvariate","lognormvariate",\r
- "expovariate","vonmisesvariate","gammavariate","triangular",\r
- "gauss","betavariate","paretovariate","weibullvariate",\r
- "getstate","setstate","jumpahead", "WichmannHill", "getrandbits",\r
- "SystemRandom"]\r
-\r
-NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)\r
-TWOPI = 2.0*_pi\r
-LOG4 = _log(4.0)\r
-SG_MAGICCONST = 1.0 + _log(4.5)\r
-BPF = 53 # Number of bits in a float\r
-RECIP_BPF = 2**-BPF\r
-\r
-\r
-# Translated by Guido van Rossum from C source provided by\r
-# Adrian Baddeley. Adapted by Raymond Hettinger for use with\r
-# the Mersenne Twister and os.urandom() core generators.\r
-\r
-import _random\r
-\r
-class Random(_random.Random):\r
- """Random number generator base class used by bound module functions.\r
-\r
- Used to instantiate instances of Random to get generators that don't\r
- share state. Especially useful for multi-threaded programs, creating\r
- a different instance of Random for each thread, and using the jumpahead()\r
- method to ensure that the generated sequences seen by each thread don't\r
- overlap.\r
-\r
- Class Random can also be subclassed if you want to use a different basic\r
- generator of your own devising: in that case, override the following\r
- methods: random(), seed(), getstate(), setstate() and jumpahead().\r
- Optionally, implement a getrandbits() method so that randrange() can cover\r
- arbitrarily large ranges.\r
-\r
- """\r
-\r
- VERSION = 3 # used by getstate/setstate\r
-\r
- def __init__(self, x=None):\r
- """Initialize an instance.\r
-\r
- Optional argument x controls seeding, as for Random.seed().\r
- """\r
-\r
- self.seed(x)\r
- self.gauss_next = None\r
-\r
- def seed(self, a=None):\r
- """Initialize internal state from hashable object.\r
-\r
- None or no argument seeds from current time or from an operating\r
- system specific randomness source if available.\r
-\r
- If a is not None or an int or long, hash(a) is used instead.\r
- """\r
-\r
- if a is None:\r
- try:\r
- # Seed with enough bytes to span the 19937 bit\r
- # state space for the Mersenne Twister\r
- a = long(_hexlify(_urandom(2500)), 16)\r
- except NotImplementedError:\r
- import time\r
- a = long(time.time() * 256) # use fractional seconds\r
-\r
- super(Random, self).seed(a)\r
- self.gauss_next = None\r
-\r
- def getstate(self):\r
- """Return internal state; can be passed to setstate() later."""\r
- return self.VERSION, super(Random, self).getstate(), self.gauss_next\r
-\r
- def setstate(self, state):\r
- """Restore internal state from object returned by getstate()."""\r
- version = state[0]\r
- if version == 3:\r
- version, internalstate, self.gauss_next = state\r
- super(Random, self).setstate(internalstate)\r
- elif version == 2:\r
- version, internalstate, self.gauss_next = state\r
- # In version 2, the state was saved as signed ints, which causes\r
- # inconsistencies between 32/64-bit systems. The state is\r
- # really unsigned 32-bit ints, so we convert negative ints from\r
- # version 2 to positive longs for version 3.\r
- try:\r
- internalstate = tuple( long(x) % (2**32) for x in internalstate )\r
- except ValueError, e:\r
- raise TypeError, e\r
- super(Random, self).setstate(internalstate)\r
- else:\r
- raise ValueError("state with version %s passed to "\r
- "Random.setstate() of version %s" %\r
- (version, self.VERSION))\r
-\r
- def jumpahead(self, n):\r
- """Change the internal state to one that is likely far away\r
- from the current state. This method will not be in Py3.x,\r
- so it is better to simply reseed.\r
- """\r
- # The super.jumpahead() method uses shuffling to change state,\r
- # so it needs a large and "interesting" n to work with. Here,\r
- # we use hashing to create a large n for the shuffle.\r
- s = repr(n) + repr(self.getstate())\r
- n = int(_hashlib.new('sha512', s).hexdigest(), 16)\r
- super(Random, self).jumpahead(n)\r
-\r
-## ---- Methods below this point do not need to be overridden when\r
-## ---- subclassing for the purpose of using a different core generator.\r
-\r
-## -------------------- pickle support -------------------\r
-\r
- def __getstate__(self): # for pickle\r
- return self.getstate()\r
-\r
- def __setstate__(self, state): # for pickle\r
- self.setstate(state)\r
-\r
- def __reduce__(self):\r
- return self.__class__, (), self.getstate()\r
-\r
-## -------------------- integer methods -------------------\r
-\r
- def randrange(self, start, stop=None, step=1, _int=int, _maxwidth=1L<<BPF):\r
- """Choose a random item from range(start, stop[, step]).\r
-\r
- This fixes the problem with randint() which includes the\r
- endpoint; in Python this is usually not what you want.\r
-\r
- """\r
-\r
- # This code is a bit messy to make it fast for the\r
- # common case while still doing adequate error checking.\r
- istart = _int(start)\r
- if istart != start:\r
- raise ValueError, "non-integer arg 1 for randrange()"\r
- if stop is None:\r
- if istart > 0:\r
- if istart >= _maxwidth:\r
- return self._randbelow(istart)\r
- return _int(self.random() * istart)\r
- raise ValueError, "empty range for randrange()"\r
-\r
- # stop argument supplied.\r
- istop = _int(stop)\r
- if istop != stop:\r
- raise ValueError, "non-integer stop for randrange()"\r
- width = istop - istart\r
- if step == 1 and width > 0:\r
- # Note that\r
- # int(istart + self.random()*width)\r
- # instead would be incorrect. For example, consider istart\r
- # = -2 and istop = 0. Then the guts would be in\r
- # -2.0 to 0.0 exclusive on both ends (ignoring that random()\r
- # might return 0.0), and because int() truncates toward 0, the\r
- # final result would be -1 or 0 (instead of -2 or -1).\r
- # istart + int(self.random()*width)\r
- # would also be incorrect, for a subtler reason: the RHS\r
- # can return a long, and then randrange() would also return\r
- # a long, but we're supposed to return an int (for backward\r
- # compatibility).\r
-\r
- if width >= _maxwidth:\r
- return _int(istart + self._randbelow(width))\r
- return _int(istart + _int(self.random()*width))\r
- if step == 1:\r
- raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width)\r
-\r
- # Non-unit step argument supplied.\r
- istep = _int(step)\r
- if istep != step:\r
- raise ValueError, "non-integer step for randrange()"\r
- if istep > 0:\r
- n = (width + istep - 1) // istep\r
- elif istep < 0:\r
- n = (width + istep + 1) // istep\r
- else:\r
- raise ValueError, "zero step for randrange()"\r
-\r
- if n <= 0:\r
- raise ValueError, "empty range for randrange()"\r
-\r
- if n >= _maxwidth:\r
- return istart + istep*self._randbelow(n)\r
- return istart + istep*_int(self.random() * n)\r
-\r
- def randint(self, a, b):\r
- """Return random integer in range [a, b], including both end points.\r
- """\r
-\r
- return self.randrange(a, b+1)\r
-\r
- def _randbelow(self, n, _log=_log, _int=int, _maxwidth=1L<<BPF,\r
- _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType):\r
- """Return a random int in the range [0,n)\r
-\r
- Handles the case where n has more bits than returned\r
- by a single call to the underlying generator.\r
- """\r
-\r
- try:\r
- getrandbits = self.getrandbits\r
- except AttributeError:\r
- pass\r
- else:\r
- # Only call self.getrandbits if the original random() builtin method\r
- # has not been overridden or if a new getrandbits() was supplied.\r
- # This assures that the two methods correspond.\r
- if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method:\r
- k = _int(1.00001 + _log(n-1, 2.0)) # 2**k > n-1 > 2**(k-2)\r
- r = getrandbits(k)\r
- while r >= n:\r
- r = getrandbits(k)\r
- return r\r
- if n >= _maxwidth:\r
- _warn("Underlying random() generator does not supply \n"\r
- "enough bits to choose from a population range this large")\r
- return _int(self.random() * n)\r
-\r
-## -------------------- sequence methods -------------------\r
-\r
- def choice(self, seq):\r
- """Choose a random element from a non-empty sequence."""\r
- return seq[int(self.random() * len(seq))] # raises IndexError if seq is empty\r
-\r
- def shuffle(self, x, random=None):\r
- """x, random=random.random -> shuffle list x in place; return None.\r
-\r
- Optional arg random is a 0-argument function returning a random\r
- float in [0.0, 1.0); by default, the standard random.random.\r
-\r
- """\r
-\r
- if random is None:\r
- random = self.random\r
- _int = int\r
- for i in reversed(xrange(1, len(x))):\r
- # pick an element in x[:i+1] with which to exchange x[i]\r
- j = _int(random() * (i+1))\r
- x[i], x[j] = x[j], x[i]\r
-\r
- def sample(self, population, k):\r
- """Chooses k unique random elements from a population sequence.\r
-\r
- Returns a new list containing elements from the population while\r
- leaving the original population unchanged. The resulting list is\r
- in selection order so that all sub-slices will also be valid random\r
- samples. This allows raffle winners (the sample) to be partitioned\r
- into grand prize and second place winners (the subslices).\r
-\r
- Members of the population need not be hashable or unique. If the\r
- population contains repeats, then each occurrence is a possible\r
- selection in the sample.\r
-\r
- To choose a sample in a range of integers, use xrange as an argument.\r
- This is especially fast and space efficient for sampling from a\r
- large population: sample(xrange(10000000), 60)\r
- """\r
-\r
- # Sampling without replacement entails tracking either potential\r
- # selections (the pool) in a list or previous selections in a set.\r
-\r
- # When the number of selections is small compared to the\r
- # population, then tracking selections is efficient, requiring\r
- # only a small set and an occasional reselection. For\r
- # a larger number of selections, the pool tracking method is\r
- # preferred since the list takes less space than the\r
- # set and it doesn't suffer from frequent reselections.\r
-\r
- n = len(population)\r
- if not 0 <= k <= n:\r
- raise ValueError("sample larger than population")\r
- random = self.random\r
- _int = int\r
- result = [None] * k\r
- setsize = 21 # size of a small set minus size of an empty list\r
- if k > 5:\r
- setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets\r
- if n <= setsize or hasattr(population, "keys"):\r
- # An n-length list is smaller than a k-length set, or this is a\r
- # mapping type so the other algorithm wouldn't work.\r
- pool = list(population)\r
- for i in xrange(k): # invariant: non-selected at [0,n-i)\r
- j = _int(random() * (n-i))\r
- result[i] = pool[j]\r
- pool[j] = pool[n-i-1] # move non-selected item into vacancy\r
- else:\r
- try:\r
- selected = set()\r
- selected_add = selected.add\r
- for i in xrange(k):\r
- j = _int(random() * n)\r
- while j in selected:\r
- j = _int(random() * n)\r
- selected_add(j)\r
- result[i] = population[j]\r
- except (TypeError, KeyError): # handle (at least) sets\r
- if isinstance(population, list):\r
- raise\r
- return self.sample(tuple(population), k)\r
- return result\r
-\r
-## -------------------- real-valued distributions -------------------\r
-\r
-## -------------------- uniform distribution -------------------\r
-\r
- def uniform(self, a, b):\r
- "Get a random number in the range [a, b) or [a, b] depending on rounding."\r
- return a + (b-a) * self.random()\r
-\r
-## -------------------- triangular --------------------\r
-\r
- def triangular(self, low=0.0, high=1.0, mode=None):\r
- """Triangular distribution.\r
-\r
- Continuous distribution bounded by given lower and upper limits,\r
- and having a given mode value in-between.\r
-\r
- http://en.wikipedia.org/wiki/Triangular_distribution\r
-\r
- """\r
- u = self.random()\r
- try:\r
- c = 0.5 if mode is None else (mode - low) / (high - low)\r
- except ZeroDivisionError:\r
- return low\r
- if u > c:\r
- u = 1.0 - u\r
- c = 1.0 - c\r
- low, high = high, low\r
- return low + (high - low) * (u * c) ** 0.5\r
-\r
-## -------------------- normal distribution --------------------\r
-\r
- def normalvariate(self, mu, sigma):\r
- """Normal distribution.\r
-\r
- mu is the mean, and sigma is the standard deviation.\r
-\r
- """\r
- # mu = mean, sigma = standard deviation\r
-\r
- # Uses Kinderman and Monahan method. Reference: Kinderman,\r
- # A.J. and Monahan, J.F., "Computer generation of random\r
- # variables using the ratio of uniform deviates", ACM Trans\r
- # Math Software, 3, (1977), pp257-260.\r
-\r
- random = self.random\r
- while 1:\r
- u1 = random()\r
- u2 = 1.0 - random()\r
- z = NV_MAGICCONST*(u1-0.5)/u2\r
- zz = z*z/4.0\r
- if zz <= -_log(u2):\r
- break\r
- return mu + z*sigma\r
-\r
-## -------------------- lognormal distribution --------------------\r
-\r
- def lognormvariate(self, mu, sigma):\r
- """Log normal distribution.\r
-\r
- If you take the natural logarithm of this distribution, you'll get a\r
- normal distribution with mean mu and standard deviation sigma.\r
- mu can have any value, and sigma must be greater than zero.\r
-\r
- """\r
- return _exp(self.normalvariate(mu, sigma))\r
-\r
-## -------------------- exponential distribution --------------------\r
-\r
- def expovariate(self, lambd):\r
- """Exponential distribution.\r
-\r
- lambd is 1.0 divided by the desired mean. It should be\r
- nonzero. (The parameter would be called "lambda", but that is\r
- a reserved word in Python.) Returned values range from 0 to\r
- positive infinity if lambd is positive, and from negative\r
- infinity to 0 if lambd is negative.\r
-\r
- """\r
- # lambd: rate lambd = 1/mean\r
- # ('lambda' is a Python reserved word)\r
-\r
- # we use 1-random() instead of random() to preclude the\r
- # possibility of taking the log of zero.\r
- return -_log(1.0 - self.random())/lambd\r
-\r
-## -------------------- von Mises distribution --------------------\r
-\r
- def vonmisesvariate(self, mu, kappa):\r
- """Circular data distribution.\r
-\r
- mu is the mean angle, expressed in radians between 0 and 2*pi, and\r
- kappa is the concentration parameter, which must be greater than or\r
- equal to zero. If kappa is equal to zero, this distribution reduces\r
- to a uniform random angle over the range 0 to 2*pi.\r
-\r
- """\r
- # mu: mean angle (in radians between 0 and 2*pi)\r
- # kappa: concentration parameter kappa (>= 0)\r
- # if kappa = 0 generate uniform random angle\r
-\r
- # Based upon an algorithm published in: Fisher, N.I.,\r
- # "Statistical Analysis of Circular Data", Cambridge\r
- # University Press, 1993.\r
-\r
- # Thanks to Magnus Kessler for a correction to the\r
- # implementation of step 4.\r
-\r
- random = self.random\r
- if kappa <= 1e-6:\r
- return TWOPI * random()\r
-\r
- s = 0.5 / kappa\r
- r = s + _sqrt(1.0 + s * s)\r
-\r
- while 1:\r
- u1 = random()\r
- z = _cos(_pi * u1)\r
-\r
- d = z / (r + z)\r
- u2 = random()\r
- if u2 < 1.0 - d * d or u2 <= (1.0 - d) * _exp(d):\r
- break\r
-\r
- q = 1.0 / r\r
- f = (q + z) / (1.0 + q * z)\r
- u3 = random()\r
- if u3 > 0.5:\r
- theta = (mu + _acos(f)) % TWOPI\r
- else:\r
- theta = (mu - _acos(f)) % TWOPI\r
-\r
- return theta\r
-\r
-## -------------------- gamma distribution --------------------\r
-\r
- def gammavariate(self, alpha, beta):\r
- """Gamma distribution. Not the gamma function!\r
-\r
- Conditions on the parameters are alpha > 0 and beta > 0.\r
-\r
- The probability distribution function is:\r
-\r
- x ** (alpha - 1) * math.exp(-x / beta)\r
- pdf(x) = --------------------------------------\r
- math.gamma(alpha) * beta ** alpha\r
-\r
- """\r
-\r
- # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2\r
-\r
- # Warning: a few older sources define the gamma distribution in terms\r
- # of alpha > -1.0\r
- if alpha <= 0.0 or beta <= 0.0:\r
- raise ValueError, 'gammavariate: alpha and beta must be > 0.0'\r
-\r
- random = self.random\r
- if alpha > 1.0:\r
-\r
- # Uses R.C.H. Cheng, "The generation of Gamma\r
- # variables with non-integral shape parameters",\r
- # Applied Statistics, (1977), 26, No. 1, p71-74\r
-\r
- ainv = _sqrt(2.0 * alpha - 1.0)\r
- bbb = alpha - LOG4\r
- ccc = alpha + ainv\r
-\r
- while 1:\r
- u1 = random()\r
- if not 1e-7 < u1 < .9999999:\r
- continue\r
- u2 = 1.0 - random()\r
- v = _log(u1/(1.0-u1))/ainv\r
- x = alpha*_exp(v)\r
- z = u1*u1*u2\r
- r = bbb+ccc*v-x\r
- if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z):\r
- return x * beta\r
-\r
- elif alpha == 1.0:\r
- # expovariate(1)\r
- u = random()\r
- while u <= 1e-7:\r
- u = random()\r
- return -_log(u) * beta\r
-\r
- else: # alpha is between 0 and 1 (exclusive)\r
-\r
- # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle\r
-\r
- while 1:\r
- u = random()\r
- b = (_e + alpha)/_e\r
- p = b*u\r
- if p <= 1.0:\r
- x = p ** (1.0/alpha)\r
- else:\r
- x = -_log((b-p)/alpha)\r
- u1 = random()\r
- if p > 1.0:\r
- if u1 <= x ** (alpha - 1.0):\r
- break\r
- elif u1 <= _exp(-x):\r
- break\r
- return x * beta\r
-\r
-## -------------------- Gauss (faster alternative) --------------------\r
-\r
- def gauss(self, mu, sigma):\r
- """Gaussian distribution.\r
-\r
- mu is the mean, and sigma is the standard deviation. This is\r
- slightly faster than the normalvariate() function.\r
-\r
- Not thread-safe without a lock around calls.\r
-\r
- """\r
-\r
- # When x and y are two variables from [0, 1), uniformly\r
- # distributed, then\r
- #\r
- # cos(2*pi*x)*sqrt(-2*log(1-y))\r
- # sin(2*pi*x)*sqrt(-2*log(1-y))\r
- #\r
- # are two *independent* variables with normal distribution\r
- # (mu = 0, sigma = 1).\r
- # (Lambert Meertens)\r
- # (corrected version; bug discovered by Mike Miller, fixed by LM)\r
-\r
- # Multithreading note: When two threads call this function\r
- # simultaneously, it is possible that they will receive the\r
- # same return value. The window is very small though. To\r
- # avoid this, you have to use a lock around all calls. (I\r
- # didn't want to slow this down in the serial case by using a\r
- # lock here.)\r
-\r
- random = self.random\r
- z = self.gauss_next\r
- self.gauss_next = None\r
- if z is None:\r
- x2pi = random() * TWOPI\r
- g2rad = _sqrt(-2.0 * _log(1.0 - random()))\r
- z = _cos(x2pi) * g2rad\r
- self.gauss_next = _sin(x2pi) * g2rad\r
-\r
- return mu + z*sigma\r
-\r
-## -------------------- beta --------------------\r
-## See\r
-## http://mail.python.org/pipermail/python-bugs-list/2001-January/003752.html\r
-## for Ivan Frohne's insightful analysis of why the original implementation:\r
-##\r
-## def betavariate(self, alpha, beta):\r
-## # Discrete Event Simulation in C, pp 87-88.\r
-##\r
-## y = self.expovariate(alpha)\r
-## z = self.expovariate(1.0/beta)\r
-## return z/(y+z)\r
-##\r
-## was dead wrong, and how it probably got that way.\r
-\r
- def betavariate(self, alpha, beta):\r
- """Beta distribution.\r
-\r
- Conditions on the parameters are alpha > 0 and beta > 0.\r
- Returned values range between 0 and 1.\r
-\r
- """\r
-\r
- # This version due to Janne Sinkkonen, and matches all the std\r
- # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").\r
- y = self.gammavariate(alpha, 1.)\r
- if y == 0:\r
- return 0.0\r
- else:\r
- return y / (y + self.gammavariate(beta, 1.))\r
-\r
-## -------------------- Pareto --------------------\r
-\r
- def paretovariate(self, alpha):\r
- """Pareto distribution. alpha is the shape parameter."""\r
- # Jain, pg. 495\r
-\r
- u = 1.0 - self.random()\r
- return 1.0 / pow(u, 1.0/alpha)\r
-\r
-## -------------------- Weibull --------------------\r
-\r
- def weibullvariate(self, alpha, beta):\r
- """Weibull distribution.\r
-\r
- alpha is the scale parameter and beta is the shape parameter.\r
-\r
- """\r
- # Jain, pg. 499; bug fix courtesy Bill Arms\r
-\r
- u = 1.0 - self.random()\r
- return alpha * pow(-_log(u), 1.0/beta)\r
-\r
-## -------------------- Wichmann-Hill -------------------\r
-\r
-class WichmannHill(Random):\r
-\r
- VERSION = 1 # used by getstate/setstate\r
-\r
- def seed(self, a=None):\r
- """Initialize internal state from hashable object.\r
-\r
- None or no argument seeds from current time or from an operating\r
- system specific randomness source if available.\r
-\r
- If a is not None or an int or long, hash(a) is used instead.\r
-\r
- If a is an int or long, a is used directly. Distinct values between\r
- 0 and 27814431486575L inclusive are guaranteed to yield distinct\r
- internal states (this guarantee is specific to the default\r
- Wichmann-Hill generator).\r
- """\r
-\r
- if a is None:\r
- try:\r
- a = long(_hexlify(_urandom(16)), 16)\r
- except NotImplementedError:\r
- import time\r
- a = long(time.time() * 256) # use fractional seconds\r
-\r
- if not isinstance(a, (int, long)):\r
- a = hash(a)\r
-\r
- a, x = divmod(a, 30268)\r
- a, y = divmod(a, 30306)\r
- a, z = divmod(a, 30322)\r
- self._seed = int(x)+1, int(y)+1, int(z)+1\r
-\r
- self.gauss_next = None\r
-\r
- def random(self):\r
- """Get the next random number in the range [0.0, 1.0)."""\r
-\r
- # Wichman-Hill random number generator.\r
- #\r
- # Wichmann, B. A. & Hill, I. D. (1982)\r
- # Algorithm AS 183:\r
- # An efficient and portable pseudo-random number generator\r
- # Applied Statistics 31 (1982) 188-190\r
- #\r
- # see also:\r
- # Correction to Algorithm AS 183\r
- # Applied Statistics 33 (1984) 123\r
- #\r
- # McLeod, A. I. (1985)\r
- # A remark on Algorithm AS 183\r
- # Applied Statistics 34 (1985),198-200\r
-\r
- # This part is thread-unsafe:\r
- # BEGIN CRITICAL SECTION\r
- x, y, z = self._seed\r
- x = (171 * x) % 30269\r
- y = (172 * y) % 30307\r
- z = (170 * z) % 30323\r
- self._seed = x, y, z\r
- # END CRITICAL SECTION\r
-\r
- # Note: on a platform using IEEE-754 double arithmetic, this can\r
- # never return 0.0 (asserted by Tim; proof too long for a comment).\r
- return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0\r
-\r
- def getstate(self):\r
- """Return internal state; can be passed to setstate() later."""\r
- return self.VERSION, self._seed, self.gauss_next\r
-\r
- def setstate(self, state):\r
- """Restore internal state from object returned by getstate()."""\r
- version = state[0]\r
- if version == 1:\r
- version, self._seed, self.gauss_next = state\r
- else:\r
- raise ValueError("state with version %s passed to "\r
- "Random.setstate() of version %s" %\r
- (version, self.VERSION))\r
-\r
- def jumpahead(self, n):\r
- """Act as if n calls to random() were made, but quickly.\r
-\r
- n is an int, greater than or equal to 0.\r
-\r
- Example use: If you have 2 threads and know that each will\r
- consume no more than a million random numbers, create two Random\r
- objects r1 and r2, then do\r
- r2.setstate(r1.getstate())\r
- r2.jumpahead(1000000)\r
- Then r1 and r2 will use guaranteed-disjoint segments of the full\r
- period.\r
- """\r
-\r
- if not n >= 0:\r
- raise ValueError("n must be >= 0")\r
- x, y, z = self._seed\r
- x = int(x * pow(171, n, 30269)) % 30269\r
- y = int(y * pow(172, n, 30307)) % 30307\r
- z = int(z * pow(170, n, 30323)) % 30323\r
- self._seed = x, y, z\r
-\r
- def __whseed(self, x=0, y=0, z=0):\r
- """Set the Wichmann-Hill seed from (x, y, z).\r
-\r
- These must be integers in the range [0, 256).\r
- """\r
-\r
- if not type(x) == type(y) == type(z) == int:\r
- raise TypeError('seeds must be integers')\r
- if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256):\r
- raise ValueError('seeds must be in range(0, 256)')\r
- if 0 == x == y == z:\r
- # Initialize from current time\r
- import time\r
- t = long(time.time() * 256)\r
- t = int((t&0xffffff) ^ (t>>24))\r
- t, x = divmod(t, 256)\r
- t, y = divmod(t, 256)\r
- t, z = divmod(t, 256)\r
- # Zero is a poor seed, so substitute 1\r
- self._seed = (x or 1, y or 1, z or 1)\r
-\r
- self.gauss_next = None\r
-\r
- def whseed(self, a=None):\r
- """Seed from hashable object's hash code.\r
-\r
- None or no argument seeds from current time. It is not guaranteed\r
- that objects with distinct hash codes lead to distinct internal\r
- states.\r
-\r
- This is obsolete, provided for compatibility with the seed routine\r
- used prior to Python 2.1. Use the .seed() method instead.\r
- """\r
-\r
- if a is None:\r
- self.__whseed()\r
- return\r
- a = hash(a)\r
- a, x = divmod(a, 256)\r
- a, y = divmod(a, 256)\r
- a, z = divmod(a, 256)\r
- x = (x + a) % 256 or 1\r
- y = (y + a) % 256 or 1\r
- z = (z + a) % 256 or 1\r
- self.__whseed(x, y, z)\r
-\r
-## --------------- Operating System Random Source ------------------\r
-\r
-class SystemRandom(Random):\r
- """Alternate random number generator using sources provided\r
- by the operating system (such as /dev/urandom on Unix or\r
- CryptGenRandom on Windows).\r
-\r
- Not available on all systems (see os.urandom() for details).\r
- """\r
-\r
- def random(self):\r
- """Get the next random number in the range [0.0, 1.0)."""\r
- return (long(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF\r
-\r
- def getrandbits(self, k):\r
- """getrandbits(k) -> x. Generates a long int with k random bits."""\r
- if k <= 0:\r
- raise ValueError('number of bits must be greater than zero')\r
- if k != int(k):\r
- raise TypeError('number of bits should be an integer')\r
- bytes = (k + 7) // 8 # bits / 8 and rounded up\r
- x = long(_hexlify(_urandom(bytes)), 16)\r
- return x >> (bytes * 8 - k) # trim excess bits\r
-\r
- def _stub(self, *args, **kwds):\r
- "Stub method. Not used for a system random number generator."\r
- return None\r
- seed = jumpahead = _stub\r
-\r
- def _notimplemented(self, *args, **kwds):\r
- "Method should not be called for a system random number generator."\r
- raise NotImplementedError('System entropy source does not have state.')\r
- getstate = setstate = _notimplemented\r
-\r
-## -------------------- test program --------------------\r
-\r
-def _test_generator(n, func, args):\r
- import time\r
- print n, 'times', func.__name__\r
- total = 0.0\r
- sqsum = 0.0\r
- smallest = 1e10\r
- largest = -1e10\r
- t0 = time.time()\r
- for i in range(n):\r
- x = func(*args)\r
- total += x\r
- sqsum = sqsum + x*x\r
- smallest = min(x, smallest)\r
- largest = max(x, largest)\r
- t1 = time.time()\r
- print round(t1-t0, 3), 'sec,',\r
- avg = total/n\r
- stddev = _sqrt(sqsum/n - avg*avg)\r
- print 'avg %g, stddev %g, min %g, max %g' % \\r
- (avg, stddev, smallest, largest)\r
-\r
-\r
-def _test(N=2000):\r
- _test_generator(N, random, ())\r
- _test_generator(N, normalvariate, (0.0, 1.0))\r
- _test_generator(N, lognormvariate, (0.0, 1.0))\r
- _test_generator(N, vonmisesvariate, (0.0, 1.0))\r
- _test_generator(N, gammavariate, (0.01, 1.0))\r
- _test_generator(N, gammavariate, (0.1, 1.0))\r
- _test_generator(N, gammavariate, (0.1, 2.0))\r
- _test_generator(N, gammavariate, (0.5, 1.0))\r
- _test_generator(N, gammavariate, (0.9, 1.0))\r
- _test_generator(N, gammavariate, (1.0, 1.0))\r
- _test_generator(N, gammavariate, (2.0, 1.0))\r
- _test_generator(N, gammavariate, (20.0, 1.0))\r
- _test_generator(N, gammavariate, (200.0, 1.0))\r
- _test_generator(N, gauss, (0.0, 1.0))\r
- _test_generator(N, betavariate, (3.0, 3.0))\r
- _test_generator(N, triangular, (0.0, 1.0, 1.0/3.0))\r
-\r
-# Create one instance, seeded from current time, and export its methods\r
-# as module-level functions. The functions share state across all uses\r
-#(both in the user's code and in the Python libraries), but that's fine\r
-# for most programs and is easier for the casual user than making them\r
-# instantiate their own Random() instance.\r
-\r
-_inst = Random()\r
-seed = _inst.seed\r
-random = _inst.random\r
-uniform = _inst.uniform\r
-triangular = _inst.triangular\r
-randint = _inst.randint\r
-choice = _inst.choice\r
-randrange = _inst.randrange\r
-sample = _inst.sample\r
-shuffle = _inst.shuffle\r
-normalvariate = _inst.normalvariate\r
-lognormvariate = _inst.lognormvariate\r
-expovariate = _inst.expovariate\r
-vonmisesvariate = _inst.vonmisesvariate\r
-gammavariate = _inst.gammavariate\r
-gauss = _inst.gauss\r
-betavariate = _inst.betavariate\r
-paretovariate = _inst.paretovariate\r
-weibullvariate = _inst.weibullvariate\r
-getstate = _inst.getstate\r
-setstate = _inst.setstate\r
-jumpahead = _inst.jumpahead\r
-getrandbits = _inst.getrandbits\r
-\r
-if __name__ == '__main__':\r
- _test()\r
projects / mirror_edk2.git / blobdiff