1"""Random variable generators. 2 3 integers 4 -------- 5 uniform within range 6 7 sequences 8 --------- 9 pick random element 10 pick random sample 11 generate random permutation 12 13 distributions on the real line: 14 ------------------------------ 15 uniform 16 triangular 17 normal (Gaussian) 18 lognormal 19 negative exponential 20 gamma 21 beta 22 pareto 23 Weibull 24 25 distributions on the circle (angles 0 to 2pi) 26 --------------------------------------------- 27 circular uniform 28 von Mises 29 30General notes on the underlying Mersenne Twister core generator: 31 32* The period is 2**19937-1. 33* It is one of the most extensively tested generators in existence. 34* Without a direct way to compute N steps forward, the semantics of 35 jumpahead(n) are weakened to simply jump to another distant state and rely 36 on the large period to avoid overlapping sequences. 37* The random() method is implemented in C, executes in a single Python step, 38 and is, therefore, threadsafe. 39 40""" 41 42from __future__ import division 43from warnings import warn as _warn 44from types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType 45from math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil 46from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin 47from os import urandom as _urandom 48from binascii import hexlify as _hexlify 49import hashlib as _hashlib 50 51__all__ = ["Random","seed","random","uniform","randint","choice","sample", 52 "randrange","shuffle","normalvariate","lognormvariate", 53 "expovariate","vonmisesvariate","gammavariate","triangular", 54 "gauss","betavariate","paretovariate","weibullvariate", 55 "getstate","setstate","jumpahead", "WichmannHill", "getrandbits", 56 "SystemRandom"] 57 58NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) 59TWOPI = 2.0*_pi 60LOG4 = _log(4.0) 61SG_MAGICCONST = 1.0 + _log(4.5) 62BPF = 53 # Number of bits in a float 63RECIP_BPF = 2**-BPF 64 65 66# Translated by Guido van Rossum from C source provided by 67# Adrian Baddeley. Adapted by Raymond Hettinger for use with 68# the Mersenne Twister and os.urandom() core generators. 69 70import _random 71 72class Random(_random.Random): 73 """Random number generator base class used by bound module functions. 74 75 Used to instantiate instances of Random to get generators that don't 76 share state. Especially useful for multi-threaded programs, creating 77 a different instance of Random for each thread, and using the jumpahead() 78 method to ensure that the generated sequences seen by each thread don't 79 overlap. 80 81 Class Random can also be subclassed if you want to use a different basic 82 generator of your own devising: in that case, override the following 83 methods: random(), seed(), getstate(), setstate() and jumpahead(). 84 Optionally, implement a getrandbits() method so that randrange() can cover 85 arbitrarily large ranges. 86 87 """ 88 89 VERSION = 3 # used by getstate/setstate 90 91 def __init__(self, x=None): 92 """Initialize an instance. 93 94 Optional argument x controls seeding, as for Random.seed(). 95 """ 96 97 self.seed(x) 98 self.gauss_next = None 99 100 def seed(self, a=None): 101 """Initialize internal state of the random number generator. 102 103 None or no argument seeds from current time or from an operating 104 system specific randomness source if available. 105 106 If a is not None or is an int or long, hash(a) is used instead. 107 Hash values for some types are nondeterministic when the 108 PYTHONHASHSEED environment variable is enabled. 109 """ 110 111 if a is None: 112 try: 113 # Seed with enough bytes to span the 19937 bit 114 # state space for the Mersenne Twister 115 a = long(_hexlify(_urandom(2500)), 16) 116 except NotImplementedError: 117 import time 118 a = long(time.time() * 256) # use fractional seconds 119 120 super(Random, self).seed(a) 121 self.gauss_next = None 122 123 def getstate(self): 124 """Return internal state; can be passed to setstate() later.""" 125 return self.VERSION, super(Random, self).getstate(), self.gauss_next 126 127 def setstate(self, state): 128 """Restore internal state from object returned by getstate().""" 129 version = state[0] 130 if version == 3: 131 version, internalstate, self.gauss_next = state 132 super(Random, self).setstate(internalstate) 133 elif version == 2: 134 version, internalstate, self.gauss_next = state 135 # In version 2, the state was saved as signed ints, which causes 136 # inconsistencies between 32/64-bit systems. The state is 137 # really unsigned 32-bit ints, so we convert negative ints from 138 # version 2 to positive longs for version 3. 139 try: 140 internalstate = tuple( long(x) % (2**32) for x in internalstate ) 141 except ValueError, e: 142 raise TypeError, e 143 super(Random, self).setstate(internalstate) 144 else: 145 raise ValueError("state with version %s passed to " 146 "Random.setstate() of version %s" % 147 (version, self.VERSION)) 148 149 def jumpahead(self, n): 150 """Change the internal state to one that is likely far away 151 from the current state. This method will not be in Py3.x, 152 so it is better to simply reseed. 153 """ 154 # The super.jumpahead() method uses shuffling to change state, 155 # so it needs a large and "interesting" n to work with. Here, 156 # we use hashing to create a large n for the shuffle. 157 s = repr(n) + repr(self.getstate()) 158 n = int(_hashlib.new('sha512', s).hexdigest(), 16) 159 super(Random, self).jumpahead(n) 160 161## ---- Methods below this point do not need to be overridden when 162## ---- subclassing for the purpose of using a different core generator. 163 164## -------------------- pickle support ------------------- 165 166 def __getstate__(self): # for pickle 167 return self.getstate() 168 169 def __setstate__(self, state): # for pickle 170 self.setstate(state) 171 172 def __reduce__(self): 173 return self.__class__, (), self.getstate() 174 175## -------------------- integer methods ------------------- 176 177 def randrange(self, start, stop=None, step=1, _int=int, _maxwidth=1L<<BPF): 178 """Choose a random item from range(start, stop[, step]). 179 180 This fixes the problem with randint() which includes the 181 endpoint; in Python this is usually not what you want. 182 183 """ 184 185 # This code is a bit messy to make it fast for the 186 # common case while still doing adequate error checking. 187 istart = _int(start) 188 if istart != start: 189 raise ValueError, "non-integer arg 1 for randrange()" 190 if stop is None: 191 if istart > 0: 192 if istart >= _maxwidth: 193 return self._randbelow(istart) 194 return _int(self.random() * istart) 195 raise ValueError, "empty range for randrange()" 196 197 # stop argument supplied. 198 istop = _int(stop) 199 if istop != stop: 200 raise ValueError, "non-integer stop for randrange()" 201 width = istop - istart 202 if step == 1 and width > 0: 203 # Note that 204 # int(istart + self.random()*width) 205 # instead would be incorrect. For example, consider istart 206 # = -2 and istop = 0. Then the guts would be in 207 # -2.0 to 0.0 exclusive on both ends (ignoring that random() 208 # might return 0.0), and because int() truncates toward 0, the 209 # final result would be -1 or 0 (instead of -2 or -1). 210 # istart + int(self.random()*width) 211 # would also be incorrect, for a subtler reason: the RHS 212 # can return a long, and then randrange() would also return 213 # a long, but we're supposed to return an int (for backward 214 # compatibility). 215 216 if width >= _maxwidth: 217 return _int(istart + self._randbelow(width)) 218 return _int(istart + _int(self.random()*width)) 219 if step == 1: 220 raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width) 221 222 # Non-unit step argument supplied. 223 istep = _int(step) 224 if istep != step: 225 raise ValueError, "non-integer step for randrange()" 226 if istep > 0: 227 n = (width + istep - 1) // istep 228 elif istep < 0: 229 n = (width + istep + 1) // istep 230 else: 231 raise ValueError, "zero step for randrange()" 232 233 if n <= 0: 234 raise ValueError, "empty range for randrange()" 235 236 if n >= _maxwidth: 237 return istart + istep*self._randbelow(n) 238 return istart + istep*_int(self.random() * n) 239 240 def randint(self, a, b): 241 """Return random integer in range [a, b], including both end points. 242 """ 243 244 return self.randrange(a, b+1) 245 246 def _randbelow(self, n, _log=_log, _int=int, _maxwidth=1L<<BPF, 247 _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType): 248 """Return a random int in the range [0,n) 249 250 Handles the case where n has more bits than returned 251 by a single call to the underlying generator. 252 """ 253 254 try: 255 getrandbits = self.getrandbits 256 except AttributeError: 257 pass 258 else: 259 # Only call self.getrandbits if the original random() builtin method 260 # has not been overridden or if a new getrandbits() was supplied. 261 # This assures that the two methods correspond. 262 if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method: 263 k = _int(1.00001 + _log(n-1, 2.0)) # 2**k > n-1 > 2**(k-2) 264 r = getrandbits(k) 265 while r >= n: 266 r = getrandbits(k) 267 return r 268 if n >= _maxwidth: 269 _warn("Underlying random() generator does not supply \n" 270 "enough bits to choose from a population range this large") 271 return _int(self.random() * n) 272 273## -------------------- sequence methods ------------------- 274 275 def choice(self, seq): 276 """Choose a random element from a non-empty sequence.""" 277 return seq[int(self.random() * len(seq))] # raises IndexError if seq is empty 278 279 def shuffle(self, x, random=None): 280 """x, random=random.random -> shuffle list x in place; return None. 281 282 Optional arg random is a 0-argument function returning a random 283 float in [0.0, 1.0); by default, the standard random.random. 284 285 """ 286 287 if random is None: 288 random = self.random 289 _int = int 290 for i in reversed(xrange(1, len(x))): 291 # pick an element in x[:i+1] with which to exchange x[i] 292 j = _int(random() * (i+1)) 293 x[i], x[j] = x[j], x[i] 294 295 def sample(self, population, k): 296 """Chooses k unique random elements from a population sequence. 297 298 Returns a new list containing elements from the population while 299 leaving the original population unchanged. The resulting list is 300 in selection order so that all sub-slices will also be valid random 301 samples. This allows raffle winners (the sample) to be partitioned 302 into grand prize and second place winners (the subslices). 303 304 Members of the population need not be hashable or unique. If the 305 population contains repeats, then each occurrence is a possible 306 selection in the sample. 307 308 To choose a sample in a range of integers, use xrange as an argument. 309 This is especially fast and space efficient for sampling from a 310 large population: sample(xrange(10000000), 60) 311 """ 312 313 # Sampling without replacement entails tracking either potential 314 # selections (the pool) in a list or previous selections in a set. 315 316 # When the number of selections is small compared to the 317 # population, then tracking selections is efficient, requiring 318 # only a small set and an occasional reselection. For 319 # a larger number of selections, the pool tracking method is 320 # preferred since the list takes less space than the 321 # set and it doesn't suffer from frequent reselections. 322 323 n = len(population) 324 if not 0 <= k <= n: 325 raise ValueError("sample larger than population") 326 random = self.random 327 _int = int 328 result = [None] * k 329 setsize = 21 # size of a small set minus size of an empty list 330 if k > 5: 331 setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets 332 if n <= setsize or hasattr(population, "keys"): 333 # An n-length list is smaller than a k-length set, or this is a 334 # mapping type so the other algorithm wouldn't work. 335 pool = list(population) 336 for i in xrange(k): # invariant: non-selected at [0,n-i) 337 j = _int(random() * (n-i)) 338 result[i] = pool[j] 339 pool[j] = pool[n-i-1] # move non-selected item into vacancy 340 else: 341 try: 342 selected = set() 343 selected_add = selected.add 344 for i in xrange(k): 345 j = _int(random() * n) 346 while j in selected: 347 j = _int(random() * n) 348 selected_add(j) 349 result[i] = population[j] 350 except (TypeError, KeyError): # handle (at least) sets 351 if isinstance(population, list): 352 raise 353 return self.sample(tuple(population), k) 354 return result 355 356## -------------------- real-valued distributions ------------------- 357 358## -------------------- uniform distribution ------------------- 359 360 def uniform(self, a, b): 361 "Get a random number in the range [a, b) or [a, b] depending on rounding." 362 return a + (b-a) * self.random() 363 364## -------------------- triangular -------------------- 365 366 def triangular(self, low=0.0, high=1.0, mode=None): 367 """Triangular distribution. 368 369 Continuous distribution bounded by given lower and upper limits, 370 and having a given mode value in-between. 371 372 http://en.wikipedia.org/wiki/Triangular_distribution 373 374 """ 375 u = self.random() 376 try: 377 c = 0.5 if mode is None else (mode - low) / (high - low) 378 except ZeroDivisionError: 379 return low 380 if u > c: 381 u = 1.0 - u 382 c = 1.0 - c 383 low, high = high, low 384 return low + (high - low) * (u * c) ** 0.5 385 386## -------------------- normal distribution -------------------- 387 388 def normalvariate(self, mu, sigma): 389 """Normal distribution. 390 391 mu is the mean, and sigma is the standard deviation. 392 393 """ 394 # mu = mean, sigma = standard deviation 395 396 # Uses Kinderman and Monahan method. Reference: Kinderman, 397 # A.J. and Monahan, J.F., "Computer generation of random 398 # variables using the ratio of uniform deviates", ACM Trans 399 # Math Software, 3, (1977), pp257-260. 400 401 random = self.random 402 while 1: 403 u1 = random() 404 u2 = 1.0 - random() 405 z = NV_MAGICCONST*(u1-0.5)/u2 406 zz = z*z/4.0 407 if zz <= -_log(u2): 408 break 409 return mu + z*sigma 410 411## -------------------- lognormal distribution -------------------- 412 413 def lognormvariate(self, mu, sigma): 414 """Log normal distribution. 415 416 If you take the natural logarithm of this distribution, you'll get a 417 normal distribution with mean mu and standard deviation sigma. 418 mu can have any value, and sigma must be greater than zero. 419 420 """ 421 return _exp(self.normalvariate(mu, sigma)) 422 423## -------------------- exponential distribution -------------------- 424 425 def expovariate(self, lambd): 426 """Exponential distribution. 427 428 lambd is 1.0 divided by the desired mean. It should be 429 nonzero. (The parameter would be called "lambda", but that is 430 a reserved word in Python.) Returned values range from 0 to 431 positive infinity if lambd is positive, and from negative 432 infinity to 0 if lambd is negative. 433 434 """ 435 # lambd: rate lambd = 1/mean 436 # ('lambda' is a Python reserved word) 437 438 # we use 1-random() instead of random() to preclude the 439 # possibility of taking the log of zero. 440 return -_log(1.0 - self.random())/lambd 441 442## -------------------- von Mises distribution -------------------- 443 444 def vonmisesvariate(self, mu, kappa): 445 """Circular data distribution. 446 447 mu is the mean angle, expressed in radians between 0 and 2*pi, and 448 kappa is the concentration parameter, which must be greater than or 449 equal to zero. If kappa is equal to zero, this distribution reduces 450 to a uniform random angle over the range 0 to 2*pi. 451 452 """ 453 # mu: mean angle (in radians between 0 and 2*pi) 454 # kappa: concentration parameter kappa (>= 0) 455 # if kappa = 0 generate uniform random angle 456 457 # Based upon an algorithm published in: Fisher, N.I., 458 # "Statistical Analysis of Circular Data", Cambridge 459 # University Press, 1993. 460 461 # Thanks to Magnus Kessler for a correction to the 462 # implementation of step 4. 463 464 random = self.random 465 if kappa <= 1e-6: 466 return TWOPI * random() 467 468 s = 0.5 / kappa 469 r = s + _sqrt(1.0 + s * s) 470 471 while 1: 472 u1 = random() 473 z = _cos(_pi * u1) 474 475 d = z / (r + z) 476 u2 = random() 477 if u2 < 1.0 - d * d or u2 <= (1.0 - d) * _exp(d): 478 break 479 480 q = 1.0 / r 481 f = (q + z) / (1.0 + q * z) 482 u3 = random() 483 if u3 > 0.5: 484 theta = (mu + _acos(f)) % TWOPI 485 else: 486 theta = (mu - _acos(f)) % TWOPI 487 488 return theta 489 490## -------------------- gamma distribution -------------------- 491 492 def gammavariate(self, alpha, beta): 493 """Gamma distribution. Not the gamma function! 494 495 Conditions on the parameters are alpha > 0 and beta > 0. 496 497 The probability distribution function is: 498 499 x ** (alpha - 1) * math.exp(-x / beta) 500 pdf(x) = -------------------------------------- 501 math.gamma(alpha) * beta ** alpha 502 503 """ 504 505 # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 506 507 # Warning: a few older sources define the gamma distribution in terms 508 # of alpha > -1.0 509 if alpha <= 0.0 or beta <= 0.0: 510 raise ValueError, 'gammavariate: alpha and beta must be > 0.0' 511 512 random = self.random 513 if alpha > 1.0: 514 515 # Uses R.C.H. Cheng, "The generation of Gamma 516 # variables with non-integral shape parameters", 517 # Applied Statistics, (1977), 26, No. 1, p71-74 518 519 ainv = _sqrt(2.0 * alpha - 1.0) 520 bbb = alpha - LOG4 521 ccc = alpha + ainv 522 523 while 1: 524 u1 = random() 525 if not 1e-7 < u1 < .9999999: 526 continue 527 u2 = 1.0 - random() 528 v = _log(u1/(1.0-u1))/ainv 529 x = alpha*_exp(v) 530 z = u1*u1*u2 531 r = bbb+ccc*v-x 532 if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 533 return x * beta 534 535 elif alpha == 1.0: 536 # expovariate(1) 537 u = random() 538 while u <= 1e-7: 539 u = random() 540 return -_log(u) * beta 541 542 else: # alpha is between 0 and 1 (exclusive) 543 544 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 545 546 while 1: 547 u = random() 548 b = (_e + alpha)/_e 549 p = b*u 550 if p <= 1.0: 551 x = p ** (1.0/alpha) 552 else: 553 x = -_log((b-p)/alpha) 554 u1 = random() 555 if p > 1.0: 556 if u1 <= x ** (alpha - 1.0): 557 break 558 elif u1 <= _exp(-x): 559 break 560 return x * beta 561 562## -------------------- Gauss (faster alternative) -------------------- 563 564 def gauss(self, mu, sigma): 565 """Gaussian distribution. 566 567 mu is the mean, and sigma is the standard deviation. This is 568 slightly faster than the normalvariate() function. 569 570 Not thread-safe without a lock around calls. 571 572 """ 573 574 # When x and y are two variables from [0, 1), uniformly 575 # distributed, then 576 # 577 # cos(2*pi*x)*sqrt(-2*log(1-y)) 578 # sin(2*pi*x)*sqrt(-2*log(1-y)) 579 # 580 # are two *independent* variables with normal distribution 581 # (mu = 0, sigma = 1). 582 # (Lambert Meertens) 583 # (corrected version; bug discovered by Mike Miller, fixed by LM) 584 585 # Multithreading note: When two threads call this function 586 # simultaneously, it is possible that they will receive the 587 # same return value. The window is very small though. To 588 # avoid this, you have to use a lock around all calls. (I 589 # didn't want to slow this down in the serial case by using a 590 # lock here.) 591 592 random = self.random 593 z = self.gauss_next 594 self.gauss_next = None 595 if z is None: 596 x2pi = random() * TWOPI 597 g2rad = _sqrt(-2.0 * _log(1.0 - random())) 598 z = _cos(x2pi) * g2rad 599 self.gauss_next = _sin(x2pi) * g2rad 600 601 return mu + z*sigma 602 603## -------------------- beta -------------------- 604## See 605## http://mail.python.org/pipermail/python-bugs-list/2001-January/003752.html 606## for Ivan Frohne's insightful analysis of why the original implementation: 607## 608## def betavariate(self, alpha, beta): 609## # Discrete Event Simulation in C, pp 87-88. 610## 611## y = self.expovariate(alpha) 612## z = self.expovariate(1.0/beta) 613## return z/(y+z) 614## 615## was dead wrong, and how it probably got that way. 616 617 def betavariate(self, alpha, beta): 618 """Beta distribution. 619 620 Conditions on the parameters are alpha > 0 and beta > 0. 621 Returned values range between 0 and 1. 622 623 """ 624 625 # This version due to Janne Sinkkonen, and matches all the std 626 # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 627 y = self.gammavariate(alpha, 1.) 628 if y == 0: 629 return 0.0 630 else: 631 return y / (y + self.gammavariate(beta, 1.)) 632 633## -------------------- Pareto -------------------- 634 635 def paretovariate(self, alpha): 636 """Pareto distribution. alpha is the shape parameter.""" 637 # Jain, pg. 495 638 639 u = 1.0 - self.random() 640 return 1.0 / pow(u, 1.0/alpha) 641 642## -------------------- Weibull -------------------- 643 644 def weibullvariate(self, alpha, beta): 645 """Weibull distribution. 646 647 alpha is the scale parameter and beta is the shape parameter. 648 649 """ 650 # Jain, pg. 499; bug fix courtesy Bill Arms 651 652 u = 1.0 - self.random() 653 return alpha * pow(-_log(u), 1.0/beta) 654 655## -------------------- Wichmann-Hill ------------------- 656 657class WichmannHill(Random): 658 659 VERSION = 1 # used by getstate/setstate 660 661 def seed(self, a=None): 662 """Initialize internal state from hashable object. 663 664 None or no argument seeds from current time or from an operating 665 system specific randomness source if available. 666 667 If a is not None or an int or long, hash(a) is used instead. 668 669 If a is an int or long, a is used directly. Distinct values between 670 0 and 27814431486575L inclusive are guaranteed to yield distinct 671 internal states (this guarantee is specific to the default 672 Wichmann-Hill generator). 673 """ 674 675 if a is None: 676 try: 677 a = long(_hexlify(_urandom(16)), 16) 678 except NotImplementedError: 679 import time 680 a = long(time.time() * 256) # use fractional seconds 681 682 if not isinstance(a, (int, long)): 683 a = hash(a) 684 685 a, x = divmod(a, 30268) 686 a, y = divmod(a, 30306) 687 a, z = divmod(a, 30322) 688 self._seed = int(x)+1, int(y)+1, int(z)+1 689 690 self.gauss_next = None 691 692 def random(self): 693 """Get the next random number in the range [0.0, 1.0).""" 694 695 # Wichman-Hill random number generator. 696 # 697 # Wichmann, B. A. & Hill, I. D. (1982) 698 # Algorithm AS 183: 699 # An efficient and portable pseudo-random number generator 700 # Applied Statistics 31 (1982) 188-190 701 # 702 # see also: 703 # Correction to Algorithm AS 183 704 # Applied Statistics 33 (1984) 123 705 # 706 # McLeod, A. I. (1985) 707 # A remark on Algorithm AS 183 708 # Applied Statistics 34 (1985),198-200 709 710 # This part is thread-unsafe: 711 # BEGIN CRITICAL SECTION 712 x, y, z = self._seed 713 x = (171 * x) % 30269 714 y = (172 * y) % 30307 715 z = (170 * z) % 30323 716 self._seed = x, y, z 717 # END CRITICAL SECTION 718 719 # Note: on a platform using IEEE-754 double arithmetic, this can 720 # never return 0.0 (asserted by Tim; proof too long for a comment). 721 return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0 722 723 def getstate(self): 724 """Return internal state; can be passed to setstate() later.""" 725 return self.VERSION, self._seed, self.gauss_next 726 727 def setstate(self, state): 728 """Restore internal state from object returned by getstate().""" 729 version = state[0] 730 if version == 1: 731 version, self._seed, self.gauss_next = state 732 else: 733 raise ValueError("state with version %s passed to " 734 "Random.setstate() of version %s" % 735 (version, self.VERSION)) 736 737 def jumpahead(self, n): 738 """Act as if n calls to random() were made, but quickly. 739 740 n is an int, greater than or equal to 0. 741 742 Example use: If you have 2 threads and know that each will 743 consume no more than a million random numbers, create two Random 744 objects r1 and r2, then do 745 r2.setstate(r1.getstate()) 746 r2.jumpahead(1000000) 747 Then r1 and r2 will use guaranteed-disjoint segments of the full 748 period. 749 """ 750 751 if not n >= 0: 752 raise ValueError("n must be >= 0") 753 x, y, z = self._seed 754 x = int(x * pow(171, n, 30269)) % 30269 755 y = int(y * pow(172, n, 30307)) % 30307 756 z = int(z * pow(170, n, 30323)) % 30323 757 self._seed = x, y, z 758 759 def __whseed(self, x=0, y=0, z=0): 760 """Set the Wichmann-Hill seed from (x, y, z). 761 762 These must be integers in the range [0, 256). 763 """ 764 765 if not type(x) == type(y) == type(z) == int: 766 raise TypeError('seeds must be integers') 767 if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): 768 raise ValueError('seeds must be in range(0, 256)') 769 if 0 == x == y == z: 770 # Initialize from current time 771 import time 772 t = long(time.time() * 256) 773 t = int((t&0xffffff) ^ (t>>24)) 774 t, x = divmod(t, 256) 775 t, y = divmod(t, 256) 776 t, z = divmod(t, 256) 777 # Zero is a poor seed, so substitute 1 778 self._seed = (x or 1, y or 1, z or 1) 779 780 self.gauss_next = None 781 782 def whseed(self, a=None): 783 """Seed from hashable object's hash code. 784 785 None or no argument seeds from current time. It is not guaranteed 786 that objects with distinct hash codes lead to distinct internal 787 states. 788 789 This is obsolete, provided for compatibility with the seed routine 790 used prior to Python 2.1. Use the .seed() method instead. 791 """ 792 793 if a is None: 794 self.__whseed() 795 return 796 a = hash(a) 797 a, x = divmod(a, 256) 798 a, y = divmod(a, 256) 799 a, z = divmod(a, 256) 800 x = (x + a) % 256 or 1 801 y = (y + a) % 256 or 1 802 z = (z + a) % 256 or 1 803 self.__whseed(x, y, z) 804 805## --------------- Operating System Random Source ------------------ 806 807class SystemRandom(Random): 808 """Alternate random number generator using sources provided 809 by the operating system (such as /dev/urandom on Unix or 810 CryptGenRandom on Windows). 811 812 Not available on all systems (see os.urandom() for details). 813 """ 814 815 def random(self): 816 """Get the next random number in the range [0.0, 1.0).""" 817 return (long(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF 818 819 def getrandbits(self, k): 820 """getrandbits(k) -> x. Generates a long int with k random bits.""" 821 if k <= 0: 822 raise ValueError('number of bits must be greater than zero') 823 if k != int(k): 824 raise TypeError('number of bits should be an integer') 825 bytes = (k + 7) // 8 # bits / 8 and rounded up 826 x = long(_hexlify(_urandom(bytes)), 16) 827 return x >> (bytes * 8 - k) # trim excess bits 828 829 def _stub(self, *args, **kwds): 830 "Stub method. Not used for a system random number generator." 831 return None 832 seed = jumpahead = _stub 833 834 def _notimplemented(self, *args, **kwds): 835 "Method should not be called for a system random number generator." 836 raise NotImplementedError('System entropy source does not have state.') 837 getstate = setstate = _notimplemented 838 839## -------------------- test program -------------------- 840 841def _test_generator(n, func, args): 842 import time 843 print n, 'times', func.__name__ 844 total = 0.0 845 sqsum = 0.0 846 smallest = 1e10 847 largest = -1e10 848 t0 = time.time() 849 for i in range(n): 850 x = func(*args) 851 total += x 852 sqsum = sqsum + x*x 853 smallest = min(x, smallest) 854 largest = max(x, largest) 855 t1 = time.time() 856 print round(t1-t0, 3), 'sec,', 857 avg = total/n 858 stddev = _sqrt(sqsum/n - avg*avg) 859 print 'avg %g, stddev %g, min %g, max %g' % \ 860 (avg, stddev, smallest, largest) 861 862 863def _test(N=2000): 864 _test_generator(N, random, ()) 865 _test_generator(N, normalvariate, (0.0, 1.0)) 866 _test_generator(N, lognormvariate, (0.0, 1.0)) 867 _test_generator(N, vonmisesvariate, (0.0, 1.0)) 868 _test_generator(N, gammavariate, (0.01, 1.0)) 869 _test_generator(N, gammavariate, (0.1, 1.0)) 870 _test_generator(N, gammavariate, (0.1, 2.0)) 871 _test_generator(N, gammavariate, (0.5, 1.0)) 872 _test_generator(N, gammavariate, (0.9, 1.0)) 873 _test_generator(N, gammavariate, (1.0, 1.0)) 874 _test_generator(N, gammavariate, (2.0, 1.0)) 875 _test_generator(N, gammavariate, (20.0, 1.0)) 876 _test_generator(N, gammavariate, (200.0, 1.0)) 877 _test_generator(N, gauss, (0.0, 1.0)) 878 _test_generator(N, betavariate, (3.0, 3.0)) 879 _test_generator(N, triangular, (0.0, 1.0, 1.0/3.0)) 880 881# Create one instance, seeded from current time, and export its methods 882# as module-level functions. The functions share state across all uses 883#(both in the user's code and in the Python libraries), but that's fine 884# for most programs and is easier for the casual user than making them 885# instantiate their own Random() instance. 886 887_inst = Random() 888seed = _inst.seed 889random = _inst.random 890uniform = _inst.uniform 891triangular = _inst.triangular 892randint = _inst.randint 893choice = _inst.choice 894randrange = _inst.randrange 895sample = _inst.sample 896shuffle = _inst.shuffle 897normalvariate = _inst.normalvariate 898lognormvariate = _inst.lognormvariate 899expovariate = _inst.expovariate 900vonmisesvariate = _inst.vonmisesvariate 901gammavariate = _inst.gammavariate 902gauss = _inst.gauss 903betavariate = _inst.betavariate 904paretovariate = _inst.paretovariate 905weibullvariate = _inst.weibullvariate 906getstate = _inst.getstate 907setstate = _inst.setstate 908jumpahead = _inst.jumpahead 909getrandbits = _inst.getrandbits 910 911if __name__ == '__main__': 912 _test() 913