0001 """Random variable generators. 0002 0003 integers 0004 -------- 0005 uniform within range 0006 0007 sequences 0008 --------- 0009 pick random element 0010 pick random sample 0011 generate random permutation 0012 0013 distributions on the real line: 0014 ------------------------------ 0015 uniform 0016 normal (Gaussian) 0017 lognormal 0018 negative exponential 0019 gamma 0020 beta 0021 pareto 0022 Weibull 0023 0024 distributions on the circle (angles 0 to 2pi) 0025 --------------------------------------------- 0026 circular uniform 0027 von Mises 0028 0029 General notes on the underlying Mersenne Twister core generator: 0030 0031 * The period is 2**19937-1. 0032 * It is one of the most extensively tested generators in existence 0033 * Without a direct way to compute N steps forward, the 0034 semantics of jumpahead(n) are weakened to simply jump 0035 to another distant state and rely on the large period 0036 to avoid overlapping sequences. 0037 * The random() method is implemented in C, executes in 0038 a single Python step, and is, therefore, threadsafe. 0039 0040 """ 0041 0042 from warnings import warn as _warn 0043 from types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType 0044 from math import log as _log, exp as _exp, pi as _pi, e as _e 0045 from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin 0046 from math import floor as _floor 0047 from os import urandom as _urandom 0048 from binascii import hexlify as _hexlify 0049 0050 __all__ = ["Random","seed","random","uniform","randint","choice","sample", 0051 "randrange","shuffle","normalvariate","lognormvariate", 0052 "expovariate","vonmisesvariate","gammavariate", 0053 "gauss","betavariate","paretovariate","weibullvariate", 0054 "getstate","setstate","jumpahead", "WichmannHill", "getrandbits", 0055 "SystemRandom"] 0056 0057 NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) 0058 TWOPI = 2.0*_pi 0059 LOG4 = _log(4.0) 0060 SG_MAGICCONST = 1.0 + _log(4.5) 0061 BPF = 53 # Number of bits in a float 0062 RECIP_BPF = 2**-BPF 0063 0064 0065 # Translated by Guido van Rossum from C source provided by 0066 # Adrian Baddeley. Adapted by Raymond Hettinger for use with 0067 # the Mersenne Twister and os.urandom() core generators. 0068 0069 import _random 0070 0071 class Random(_random.Random): 0072 """Random number generator base class used by bound module functions. 0073 0074 Used to instantiate instances of Random to get generators that don't 0075 share state. Especially useful for multi-threaded programs, creating 0076 a different instance of Random for each thread, and using the jumpahead() 0077 method to ensure that the generated sequences seen by each thread don't 0078 overlap. 0079 0080 Class Random can also be subclassed if you want to use a different basic 0081 generator of your own devising: in that case, override the following 0082 methods: random(), seed(), getstate(), setstate() and jumpahead(). 0083 Optionally, implement a getrandombits() method so that randrange() 0084 can cover arbitrarily large ranges. 0085 0086 """ 0087 0088 VERSION = 2 # used by getstate/setstate 0089 0090 def __init__(self, x=None): 0091 """Initialize an instance. 0092 0093 Optional argument x controls seeding, as for Random.seed(). 0094 """ 0095 0096 self.seed(x) 0097 self.gauss_next = None 0098 0099 def seed(self, a=None): 0100 """Initialize internal state from hashable object. 0101 0102 None or no argument seeds from current time or from an operating 0103 system specific randomness source if available. 0104 0105 If a is not None or an int or long, hash(a) is used instead. 0106 """ 0107 0108 if a is None: 0109 try: 0110 a = long(_hexlify(_urandom(16)), 16) 0111 except NotImplementedError: 0112 import time 0113 a = long(time.time() * 256) # use fractional seconds 0114 0115 super(Random, self).seed(a) 0116 self.gauss_next = None 0117 0118 def getstate(self): 0119 """Return internal state; can be passed to setstate() later.""" 0120 return self.VERSION, super(Random, self).getstate(), self.gauss_next 0121 0122 def setstate(self, state): 0123 """Restore internal state from object returned by getstate().""" 0124 version = state[0] 0125 if version == 2: 0126 version, internalstate, self.gauss_next = state 0127 super(Random, self).setstate(internalstate) 0128 else: 0129 raise ValueError("state with version %s passed to " 0130 "Random.setstate() of version %s" % 0131 (version, self.VERSION)) 0132 0133 ## ---- Methods below this point do not need to be overridden when 0134 ## ---- subclassing for the purpose of using a different core generator. 0135 0136 ## -------------------- pickle support ------------------- 0137 0138 def __getstate__(self): # for pickle 0139 return self.getstate() 0140 0141 def __setstate__(self, state): # for pickle 0142 self.setstate(state) 0143 0144 def __reduce__(self): 0145 return self.__class__, (), self.getstate() 0146 0147 ## -------------------- integer methods ------------------- 0148 0149 def randrange(self, start, stop=None, step=1, int=int, default=None, 0150 maxwidth=1L<<BPF): 0151 """Choose a random item from range(start, stop[, step]). 0152 0153 This fixes the problem with randint() which includes the 0154 endpoint; in Python this is usually not what you want. 0155 Do not supply the 'int', 'default', and 'maxwidth' arguments. 0156 """ 0157 0158 # This code is a bit messy to make it fast for the 0159 # common case while still doing adequate error checking. 0160 istart = int(start) 0161 if istart != start: 0162 raise ValueError, "non-integer arg 1 for randrange()" 0163 if stop is default: 0164 if istart > 0: 0165 if istart >= maxwidth: 0166 return self._randbelow(istart) 0167 return int(self.random() * istart) 0168 raise ValueError, "empty range for randrange()" 0169 0170 # stop argument supplied. 0171 istop = int(stop) 0172 if istop != stop: 0173 raise ValueError, "non-integer stop for randrange()" 0174 width = istop - istart 0175 if step == 1 and width > 0: 0176 # Note that 0177 # int(istart + self.random()*width) 0178 # instead would be incorrect. For example, consider istart 0179 # = -2 and istop = 0. Then the guts would be in 0180 # -2.0 to 0.0 exclusive on both ends (ignoring that random() 0181 # might return 0.0), and because int() truncates toward 0, the 0182 # final result would be -1 or 0 (instead of -2 or -1). 0183 # istart + int(self.random()*width) 0184 # would also be incorrect, for a subtler reason: the RHS 0185 # can return a long, and then randrange() would also return 0186 # a long, but we're supposed to return an int (for backward 0187 # compatibility). 0188 0189 if width >= maxwidth: 0190 return int(istart + self._randbelow(width)) 0191 return int(istart + int(self.random()*width)) 0192 if step == 1: 0193 raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width) 0194 0195 # Non-unit step argument supplied. 0196 istep = int(step) 0197 if istep != step: 0198 raise ValueError, "non-integer step for randrange()" 0199 if istep > 0: 0200 n = (width + istep - 1) // istep 0201 elif istep < 0: 0202 n = (width + istep + 1) // istep 0203 else: 0204 raise ValueError, "zero step for randrange()" 0205 0206 if n <= 0: 0207 raise ValueError, "empty range for randrange()" 0208 0209 if n >= maxwidth: 0210 return istart + self._randbelow(n) 0211 return istart + istep*int(self.random() * n) 0212 0213 def randint(self, a, b): 0214 """Return random integer in range [a, b], including both end points. 0215 """ 0216 0217 return self.randrange(a, b+1) 0218 0219 def _randbelow(self, n, _log=_log, int=int, _maxwidth=1L<<BPF, 0220 _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType): 0221 """Return a random int in the range [0,n) 0222 0223 Handles the case where n has more bits than returned 0224 by a single call to the underlying generator. 0225 """ 0226 0227 try: 0228 getrandbits = self.getrandbits 0229 except AttributeError: 0230 pass 0231 else: 0232 # Only call self.getrandbits if the original random() builtin method 0233 # has not been overridden or if a new getrandbits() was supplied. 0234 # This assures that the two methods correspond. 0235 if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method: 0236 k = int(1.00001 + _log(n-1, 2.0)) # 2**k > n-1 > 2**(k-2) 0237 r = getrandbits(k) 0238 while r >= n: 0239 r = getrandbits(k) 0240 return r 0241 if n >= _maxwidth: 0242 _warn("Underlying random() generator does not supply \n" 0243 "enough bits to choose from a population range this large") 0244 return int(self.random() * n) 0245 0246 ## -------------------- sequence methods ------------------- 0247 0248 def choice(self, seq): 0249 """Choose a random element from a non-empty sequence.""" 0250 return seq[int(self.random() * len(seq))] # raises IndexError if seq is empty 0251 0252 def shuffle(self, x, random=None, int=int): 0253 """x, random=random.random -> shuffle list x in place; return None. 0254 0255 Optional arg random is a 0-argument function returning a random 0256 float in [0.0, 1.0); by default, the standard random.random. 0257 0258 Note that for even rather small len(x), the total number of 0259 permutations of x is larger than the period of most random number 0260 generators; this implies that "most" permutations of a long 0261 sequence can never be generated. 0262 """ 0263 0264 if random is None: 0265 random = self.random 0266 for i in reversed(xrange(1, len(x))): 0267 # pick an element in x[:i+1] with which to exchange x[i] 0268 j = int(random() * (i+1)) 0269 x[i], x[j] = x[j], x[i] 0270 0271 def sample(self, population, k): 0272 """Chooses k unique random elements from a population sequence. 0273 0274 Returns a new list containing elements from the population while 0275 leaving the original population unchanged. The resulting list is 0276 in selection order so that all sub-slices will also be valid random 0277 samples. This allows raffle winners (the sample) to be partitioned 0278 into grand prize and second place winners (the subslices). 0279 0280 Members of the population need not be hashable or unique. If the 0281 population contains repeats, then each occurrence is a possible 0282 selection in the sample. 0283 0284 To choose a sample in a range of integers, use xrange as an argument. 0285 This is especially fast and space efficient for sampling from a 0286 large population: sample(xrange(10000000), 60) 0287 """ 0288 0289 # Sampling without replacement entails tracking either potential 0290 # selections (the pool) in a list or previous selections in a 0291 # dictionary. 0292 0293 # When the number of selections is small compared to the 0294 # population, then tracking selections is efficient, requiring 0295 # only a small dictionary and an occasional reselection. For 0296 # a larger number of selections, the pool tracking method is 0297 # preferred since the list takes less space than the 0298 # dictionary and it doesn't suffer from frequent reselections. 0299 0300 n = len(population) 0301 if not 0 <= k <= n: 0302 raise ValueError, "sample larger than population" 0303 random = self.random 0304 _int = int 0305 result = [None] * k 0306 if n < 6 * k: # if n len list takes less space than a k len dict 0307 pool = list(population) 0308 for i in xrange(k): # invariant: non-selected at [0,n-i) 0309 j = _int(random() * (n-i)) 0310 result[i] = pool[j] 0311 pool[j] = pool[n-i-1] # move non-selected item into vacancy 0312 else: 0313 try: 0314 n > 0 and (population[0], population[n//2], population[n-1]) 0315 except (TypeError, KeyError): # handle sets and dictionaries 0316 population = tuple(population) 0317 selected = {} 0318 for i in xrange(k): 0319 j = _int(random() * n) 0320 while j in selected: 0321 j = _int(random() * n) 0322 result[i] = selected[j] = population[j] 0323 return result 0324 0325 ## -------------------- real-valued distributions ------------------- 0326 0327 ## -------------------- uniform distribution ------------------- 0328 0329 def uniform(self, a, b): 0330 """Get a random number in the range [a, b).""" 0331 return a + (b-a) * self.random() 0332 0333 ## -------------------- normal distribution -------------------- 0334 0335 def normalvariate(self, mu, sigma): 0336 """Normal distribution. 0337 0338 mu is the mean, and sigma is the standard deviation. 0339 0340 """ 0341 # mu = mean, sigma = standard deviation 0342 0343 # Uses Kinderman and Monahan method. Reference: Kinderman, 0344 # A.J. and Monahan, J.F., "Computer generation of random 0345 # variables using the ratio of uniform deviates", ACM Trans 0346 # Math Software, 3, (1977), pp257-260. 0347 0348 random = self.random 0349 while True: 0350 u1 = random() 0351 u2 = 1.0 - random() 0352 z = NV_MAGICCONST*(u1-0.5)/u2 0353 zz = z*z/4.0 0354 if zz <= -_log(u2): 0355 break 0356 return mu + z*sigma 0357 0358 ## -------------------- lognormal distribution -------------------- 0359 0360 def lognormvariate(self, mu, sigma): 0361 """Log normal distribution. 0362 0363 If you take the natural logarithm of this distribution, you'll get a 0364 normal distribution with mean mu and standard deviation sigma. 0365 mu can have any value, and sigma must be greater than zero. 0366 0367 """ 0368 return _exp(self.normalvariate(mu, sigma)) 0369 0370 ## -------------------- exponential distribution -------------------- 0371 0372 def expovariate(self, lambd): 0373 """Exponential distribution. 0374 0375 lambd is 1.0 divided by the desired mean. (The parameter would be 0376 called "lambda", but that is a reserved word in Python.) Returned 0377 values range from 0 to positive infinity. 0378 0379 """ 0380 # lambd: rate lambd = 1/mean 0381 # ('lambda' is a Python reserved word) 0382 0383 random = self.random 0384 u = random() 0385 while u <= 1e-7: 0386 u = random() 0387 return -_log(u)/lambd 0388 0389 ## -------------------- von Mises distribution -------------------- 0390 0391 def vonmisesvariate(self, mu, kappa): 0392 """Circular data distribution. 0393 0394 mu is the mean angle, expressed in radians between 0 and 2*pi, and 0395 kappa is the concentration parameter, which must be greater than or 0396 equal to zero. If kappa is equal to zero, this distribution reduces 0397 to a uniform random angle over the range 0 to 2*pi. 0398 0399 """ 0400 # mu: mean angle (in radians between 0 and 2*pi) 0401 # kappa: concentration parameter kappa (>= 0) 0402 # if kappa = 0 generate uniform random angle 0403 0404 # Based upon an algorithm published in: Fisher, N.I., 0405 # "Statistical Analysis of Circular Data", Cambridge 0406 # University Press, 1993. 0407 0408 # Thanks to Magnus Kessler for a correction to the 0409 # implementation of step 4. 0410 0411 random = self.random 0412 if kappa <= 1e-6: 0413 return TWOPI * random() 0414 0415 a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa) 0416 b = (a - _sqrt(2.0 * a))/(2.0 * kappa) 0417 r = (1.0 + b * b)/(2.0 * b) 0418 0419 while True: 0420 u1 = random() 0421 0422 z = _cos(_pi * u1) 0423 f = (1.0 + r * z)/(r + z) 0424 c = kappa * (r - f) 0425 0426 u2 = random() 0427 0428 if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)): 0429 break 0430 0431 u3 = random() 0432 if u3 > 0.5: 0433 theta = (mu % TWOPI) + _acos(f) 0434 else: 0435 theta = (mu % TWOPI) - _acos(f) 0436 0437 return theta 0438 0439 ## -------------------- gamma distribution -------------------- 0440 0441 def gammavariate(self, alpha, beta): 0442 """Gamma distribution. Not the gamma function! 0443 0444 Conditions on the parameters are alpha > 0 and beta > 0. 0445 0446 """ 0447 0448 # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 0449 0450 # Warning: a few older sources define the gamma distribution in terms 0451 # of alpha > -1.0 0452 if alpha <= 0.0 or beta <= 0.0: 0453 raise ValueError, 'gammavariate: alpha and beta must be > 0.0' 0454 0455 random = self.random 0456 if alpha > 1.0: 0457 0458 # Uses R.C.H. Cheng, "The generation of Gamma 0459 # variables with non-integral shape parameters", 0460 # Applied Statistics, (1977), 26, No. 1, p71-74 0461 0462 ainv = _sqrt(2.0 * alpha - 1.0) 0463 bbb = alpha - LOG4 0464 ccc = alpha + ainv 0465 0466 while True: 0467 u1 = random() 0468 if not 1e-7 < u1 < .9999999: 0469 continue 0470 u2 = 1.0 - random() 0471 v = _log(u1/(1.0-u1))/ainv 0472 x = alpha*_exp(v) 0473 z = u1*u1*u2 0474 r = bbb+ccc*v-x 0475 if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 0476 return x * beta 0477 0478 elif alpha == 1.0: 0479 # expovariate(1) 0480 u = random() 0481 while u <= 1e-7: 0482 u = random() 0483 return -_log(u) * beta 0484 0485 else: # alpha is between 0 and 1 (exclusive) 0486 0487 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 0488 0489 while True: 0490 u = random() 0491 b = (_e + alpha)/_e 0492 p = b*u 0493 if p <= 1.0: 0494 x = pow(p, 1.0/alpha) 0495 else: 0496 # p > 1 0497 x = -_log((b-p)/alpha) 0498 u1 = random() 0499 if not (((p <= 1.0) and (u1 > _exp(-x))) or 0500 ((p > 1) and (u1 > pow(x, alpha - 1.0)))): 0501 break 0502 return x * beta 0503 0504 ## -------------------- Gauss (faster alternative) -------------------- 0505 0506 def gauss(self, mu, sigma): 0507 """Gaussian distribution. 0508 0509 mu is the mean, and sigma is the standard deviation. This is 0510 slightly faster than the normalvariate() function. 0511 0512 Not thread-safe without a lock around calls. 0513 0514 """ 0515 0516 # When x and y are two variables from [0, 1), uniformly 0517 # distributed, then 0518 # 0519 # cos(2*pi*x)*sqrt(-2*log(1-y)) 0520 # sin(2*pi*x)*sqrt(-2*log(1-y)) 0521 # 0522 # are two *independent* variables with normal distribution 0523 # (mu = 0, sigma = 1). 0524 # (Lambert Meertens) 0525 # (corrected version; bug discovered by Mike Miller, fixed by LM) 0526 0527 # Multithreading note: When two threads call this function 0528 # simultaneously, it is possible that they will receive the 0529 # same return value. The window is very small though. To 0530 # avoid this, you have to use a lock around all calls. (I 0531 # didn't want to slow this down in the serial case by using a 0532 # lock here.) 0533 0534 random = self.random 0535 z = self.gauss_next 0536 self.gauss_next = None 0537 if z is None: 0538 x2pi = random() * TWOPI 0539 g2rad = _sqrt(-2.0 * _log(1.0 - random())) 0540 z = _cos(x2pi) * g2rad 0541 self.gauss_next = _sin(x2pi) * g2rad 0542 0543 return mu + z*sigma 0544 0545 ## -------------------- beta -------------------- 0546 ## See 0547 ## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470 0548 ## for Ivan Frohne's insightful analysis of why the original implementation: 0549 ## 0550 ## def betavariate(self, alpha, beta): 0551 ## # Discrete Event Simulation in C, pp 87-88. 0552 ## 0553 ## y = self.expovariate(alpha) 0554 ## z = self.expovariate(1.0/beta) 0555 ## return z/(y+z) 0556 ## 0557 ## was dead wrong, and how it probably got that way. 0558 0559 def betavariate(self, alpha, beta): 0560 """Beta distribution. 0561 0562 Conditions on the parameters are alpha > -1 and beta} > -1. 0563 Returned values range between 0 and 1. 0564 0565 """ 0566 0567 # This version due to Janne Sinkkonen, and matches all the std 0568 # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 0569 y = self.gammavariate(alpha, 1.) 0570 if y == 0: 0571 return 0.0 0572 else: 0573 return y / (y + self.gammavariate(beta, 1.)) 0574 0575 ## -------------------- Pareto -------------------- 0576 0577 def paretovariate(self, alpha): 0578 """Pareto distribution. alpha is the shape parameter.""" 0579 # Jain, pg. 495 0580 0581 u = 1.0 - self.random() 0582 return 1.0 / pow(u, 1.0/alpha) 0583 0584 ## -------------------- Weibull -------------------- 0585 0586 def weibullvariate(self, alpha, beta): 0587 """Weibull distribution. 0588 0589 alpha is the scale parameter and beta is the shape parameter. 0590 0591 """ 0592 # Jain, pg. 499; bug fix courtesy Bill Arms 0593 0594 u = 1.0 - self.random() 0595 return alpha * pow(-_log(u), 1.0/beta) 0596 0597 ## -------------------- Wichmann-Hill ------------------- 0598 0599 class WichmannHill(Random): 0600 0601 VERSION = 1 # used by getstate/setstate 0602 0603 def seed(self, a=None): 0604 """Initialize internal state from hashable object. 0605 0606 None or no argument seeds from current time or from an operating 0607 system specific randomness source if available. 0608 0609 If a is not None or an int or long, hash(a) is used instead. 0610 0611 If a is an int or long, a is used directly. Distinct values between 0612 0 and 27814431486575L inclusive are guaranteed to yield distinct 0613 internal states (this guarantee is specific to the default 0614 Wichmann-Hill generator). 0615 """ 0616 0617 if a is None: 0618 try: 0619 a = long(_hexlify(_urandom(16)), 16) 0620 except NotImplementedError: 0621 import time 0622 a = long(time.time() * 256) # use fractional seconds 0623 0624 if not isinstance(a, (int, long)): 0625 a = hash(a) 0626 0627 a, x = divmod(a, 30268) 0628 a, y = divmod(a, 30306) 0629 a, z = divmod(a, 30322) 0630 self._seed = int(x)+1, int(y)+1, int(z)+1 0631 0632 self.gauss_next = None 0633 0634 def random(self): 0635 """Get the next random number in the range [0.0, 1.0).""" 0636 0637 # Wichman-Hill random number generator. 0638 # 0639 # Wichmann, B. A. & Hill, I. D. (1982) 0640 # Algorithm AS 183: 0641 # An efficient and portable pseudo-random number generator 0642 # Applied Statistics 31 (1982) 188-190 0643 # 0644 # see also: 0645 # Correction to Algorithm AS 183 0646 # Applied Statistics 33 (1984) 123 0647 # 0648 # McLeod, A. I. (1985) 0649 # A remark on Algorithm AS 183 0650 # Applied Statistics 34 (1985),198-200 0651 0652 # This part is thread-unsafe: 0653 # BEGIN CRITICAL SECTION 0654 x, y, z = self._seed 0655 x = (171 * x) % 30269 0656 y = (172 * y) % 30307 0657 z = (170 * z) % 30323 0658 self._seed = x, y, z 0659 # END CRITICAL SECTION 0660 0661 # Note: on a platform using IEEE-754 double arithmetic, this can 0662 # never return 0.0 (asserted by Tim; proof too long for a comment). 0663 return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0 0664 0665 def getstate(self): 0666 """Return internal state; can be passed to setstate() later.""" 0667 return self.VERSION, self._seed, self.gauss_next 0668 0669 def setstate(self, state): 0670 """Restore internal state from object returned by getstate().""" 0671 version = state[0] 0672 if version == 1: 0673 version, self._seed, self.gauss_next = state 0674 else: 0675 raise ValueError("state with version %s passed to " 0676 "Random.setstate() of version %s" % 0677 (version, self.VERSION)) 0678 0679 def jumpahead(self, n): 0680 """Act as if n calls to random() were made, but quickly. 0681 0682 n is an int, greater than or equal to 0. 0683 0684 Example use: If you have 2 threads and know that each will 0685 consume no more than a million random numbers, create two Random 0686 objects r1 and r2, then do 0687 r2.setstate(r1.getstate()) 0688 r2.jumpahead(1000000) 0689 Then r1 and r2 will use guaranteed-disjoint segments of the full 0690 period. 0691 """ 0692 0693 if not n >= 0: 0694 raise ValueError("n must be >= 0") 0695 x, y, z = self._seed 0696 x = int(x * pow(171, n, 30269)) % 30269 0697 y = int(y * pow(172, n, 30307)) % 30307 0698 z = int(z * pow(170, n, 30323)) % 30323 0699 self._seed = x, y, z 0700 0701 def __whseed(self, x=0, y=0, z=0): 0702 """Set the Wichmann-Hill seed from (x, y, z). 0703 0704 These must be integers in the range [0, 256). 0705 """ 0706 0707 if not type(x) == type(y) == type(z) == int: 0708 raise TypeError('seeds must be integers') 0709 if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): 0710 raise ValueError('seeds must be in range(0, 256)') 0711 if 0 == x == y == z: 0712 # Initialize from current time 0713 import time 0714 t = long(time.time() * 256) 0715 t = int((t&0xffffff) ^ (t>>24)) 0716 t, x = divmod(t, 256) 0717 t, y = divmod(t, 256) 0718 t, z = divmod(t, 256) 0719 # Zero is a poor seed, so substitute 1 0720 self._seed = (x or 1, y or 1, z or 1) 0721 0722 self.gauss_next = None 0723 0724 def whseed(self, a=None): 0725 """Seed from hashable object's hash code. 0726 0727 None or no argument seeds from current time. It is not guaranteed 0728 that objects with distinct hash codes lead to distinct internal 0729 states. 0730 0731 This is obsolete, provided for compatibility with the seed routine 0732 used prior to Python 2.1. Use the .seed() method instead. 0733 """ 0734 0735 if a is None: 0736 self.__whseed() 0737 return 0738 a = hash(a) 0739 a, x = divmod(a, 256) 0740 a, y = divmod(a, 256) 0741 a, z = divmod(a, 256) 0742 x = (x + a) % 256 or 1 0743 y = (y + a) % 256 or 1 0744 z = (z + a) % 256 or 1 0745 self.__whseed(x, y, z) 0746 0747 ## --------------- Operating System Random Source ------------------ 0748 0749 class SystemRandom(Random): 0750 """Alternate random number generator using sources provided 0751 by the operating system (such as /dev/urandom on Unix or 0752 CryptGenRandom on Windows). 0753 0754 Not available on all systems (see os.urandom() for details). 0755 """ 0756 0757 def random(self): 0758 """Get the next random number in the range [0.0, 1.0).""" 0759 return (long(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF 0760 0761 def getrandbits(self, k): 0762 """getrandbits(k) -> x. Generates a long int with k random bits.""" 0763 if k <= 0: 0764 raise ValueError('number of bits must be greater than zero') 0765 if k != int(k): 0766 raise TypeError('number of bits should be an integer') 0767 bytes = (k + 7) // 8 # bits / 8 and rounded up 0768 x = long(_hexlify(_urandom(bytes)), 16) 0769 return x >> (bytes * 8 - k) # trim excess bits 0770 0771 def _stub(self, *args, **kwds): 0772 "Stub method. Not used for a system random number generator." 0773 return None 0774 seed = jumpahead = _stub 0775 0776 def _notimplemented(self, *args, **kwds): 0777 "Method should not be called for a system random number generator." 0778 raise NotImplementedError('System entropy source does not have state.') 0779 getstate = setstate = _notimplemented 0780 0781 ## -------------------- test program -------------------- 0782 0783 def _test_generator(n, func, args): 0784 import time 0785 print n, 'times', func.__name__ 0786 total = 0.0 0787 sqsum = 0.0 0788 smallest = 1e10 0789 largest = -1e10 0790 t0 = time.time() 0791 for i in range(n): 0792 x = func(*args) 0793 total += x 0794 sqsum = sqsum + x*x 0795 smallest = min(x, smallest) 0796 largest = max(x, largest) 0797 t1 = time.time() 0798 print round(t1-t0, 3), 'sec,', 0799 avg = total/n 0800 stddev = _sqrt(sqsum/n - avg*avg) 0801 print 'avg %g, stddev %g, min %g, max %g' % \ 0802 (avg, stddev, smallest, largest) 0803 0804 0805 def _test(N=2000): 0806 _test_generator(N, random, ()) 0807 _test_generator(N, normalvariate, (0.0, 1.0)) 0808 _test_generator(N, lognormvariate, (0.0, 1.0)) 0809 _test_generator(N, vonmisesvariate, (0.0, 1.0)) 0810 _test_generator(N, gammavariate, (0.01, 1.0)) 0811 _test_generator(N, gammavariate, (0.1, 1.0)) 0812 _test_generator(N, gammavariate, (0.1, 2.0)) 0813 _test_generator(N, gammavariate, (0.5, 1.0)) 0814 _test_generator(N, gammavariate, (0.9, 1.0)) 0815 _test_generator(N, gammavariate, (1.0, 1.0)) 0816 _test_generator(N, gammavariate, (2.0, 1.0)) 0817 _test_generator(N, gammavariate, (20.0, 1.0)) 0818 _test_generator(N, gammavariate, (200.0, 1.0)) 0819 _test_generator(N, gauss, (0.0, 1.0)) 0820 _test_generator(N, betavariate, (3.0, 3.0)) 0821 0822 # Create one instance, seeded from current time, and export its methods 0823 # as module-level functions. The functions share state across all uses 0824 #(both in the user's code and in the Python libraries), but that's fine 0825 # for most programs and is easier for the casual user than making them 0826 # instantiate their own Random() instance. 0827 0828 _inst = Random() 0829 seed = _inst.seed 0830 random = _inst.random 0831 uniform = _inst.uniform 0832 randint = _inst.randint 0833 choice = _inst.choice 0834 randrange = _inst.randrange 0835 sample = _inst.sample 0836 shuffle = _inst.shuffle 0837 normalvariate = _inst.normalvariate 0838 lognormvariate = _inst.lognormvariate 0839 expovariate = _inst.expovariate 0840 vonmisesvariate = _inst.vonmisesvariate 0841 gammavariate = _inst.gammavariate 0842 gauss = _inst.gauss 0843 betavariate = _inst.betavariate 0844 paretovariate = _inst.paretovariate 0845 weibullvariate = _inst.weibullvariate 0846 getstate = _inst.getstate 0847 setstate = _inst.setstate 0848 jumpahead = _inst.jumpahead 0849 getrandbits = _inst.getrandbits 0850 0851 if __name__ == '__main__': 0852 _test() 0853
Generated by PyXR 0.9.4