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