random

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 normal (Gaussian) 17 lognormal 18 negative exponential 19 gamma 20 beta 21 pareto 22 Weibull 23 24 distributions on the circle (angles 0 to 2pi) 25 --------------------------------------------- 26 circular uniform 27 von Mises 28 29 General notes on the underlying Mersenne Twister core generator: 30 31 * The period is 2**19937-1. 32 * It is one of the most extensively tested generators in existence. 33 * Without a direct way to compute N steps forward, the semantics of 34 jumpahead(n) are weakened to simply jump to another distant state and rely 35 on the large period to avoid overlapping sequences. 36 * The random() method is implemented in C, executes in a single Python step, 37 and is, therefore, threadsafe. 38 39 """ 40 41 from warnings import warn as _warn 42 from types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType 43 from math import log as _log, exp as _exp, pi as _pi, e as _e 44 from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin 45 from os import urandom as _urandom 46 from binascii import hexlify as _hexlify 47 48 __all__ = ["Random","seed","random","uniform","randint","choice","sample", 49 "randrange","shuffle","normalvariate","lognormvariate", 50 "expovariate","vonmisesvariate","gammavariate", 51 "gauss","betavariate","paretovariate","weibullvariate", 52 "getstate","setstate","jumpahead", "WichmannHill", "getrandbits", 53 "SystemRandom"] 54 55 NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) 56 TWOPI = 2.0*_pi 57 LOG4 = _log(4.0) 58 SG_MAGICCONST = 1.0 + _log(4.5) 59 BPF = 53 # Number of bits in a float 60 RECIP_BPF = 2**-BPF 61 62 63 # Translated by Guido van Rossum from C source provided by 64 # Adrian Baddeley. Adapted by Raymond Hettinger for use with 65 # the Mersenne Twister and os.urandom() core generators. 66 67 import _random 68

69 -class Random(_random.Random):

70 """Random number generator base class used by bound module functions. 71 72 Used to instantiate instances of Random to get generators that don't 73 share state. Especially useful for multi-threaded programs, creating 74 a different instance of Random for each thread, and using the jumpahead() 75 method to ensure that the generated sequences seen by each thread don't 76 overlap. 77 78 Class Random can also be subclassed if you want to use a different basic 79 generator of your own devising: in that case, override the following 80 methods: random(), seed(), getstate(), setstate() and jumpahead(). 81 Optionally, implement a getrandombits() method so that randrange() 82 can cover arbitrarily large ranges. 83 84 """ 85 86 VERSION = 2 # used by getstate/setstate 87

88 - def __init__(self, x=None):

89 """Initialize an instance. 90 91 Optional argument x controls seeding, as for Random.seed(). 92 """ 93 94 self.seed(x) 95 self.gauss_next = None

96

97 - def seed(self, a=None):

98 """Initialize internal state from hashable object. 99 100 None or no argument seeds from current time or from an operating 101 system specific randomness source if available. 102 103 If a is not None or an int or long, hash(a) is used instead. 104 """ 105 106 if a is None: 107 try: 108 a = long(_hexlify(_urandom(16)), 16) 109 except NotImplementedError: 110 import time 111 a = long(time.time() * 256) # use fractional seconds 112 113 super(Random, self).seed(a) 114 self.gauss_next = None

115

116 - def getstate(self):

117 """Return internal state; can be passed to setstate() later.""" 118 return self.VERSION, super(Random, self).getstate(), self.gauss_next

119

120 - def setstate(self, state):

121 """Restore internal state from object returned by getstate().""" 122 version = state[0] 123 if version == 2: 124 version, internalstate, self.gauss_next = state 125 super(Random, self).setstate(internalstate) 126 else: 127 raise ValueError("state with version %s passed to " 128 "Random.setstate() of version %s" % 129 (version, self.VERSION))

130 131 ## ---- Methods below this point do not need to be overridden when 132 ## ---- subclassing for the purpose of using a different core generator. 133 134 ## -------------------- pickle support ------------------- 135

136 - def __getstate__(self): # for pickle

137 return self.getstate()

138

139 - def __setstate__(self, state): # for pickle

140 self.setstate(state) 141

142 - def __reduce__(self):

143 return self.__class__, (), self.getstate()

144 145 ## -------------------- integer methods ------------------- 146

147 - def randrange(self, start, stop=None, step=1, int=int, default=None, 148 maxwidth=1L<<BPF):

149 """Choose a random item from range(start, stop[, step]). 150 151 This fixes the problem with randint() which includes the 152 endpoint; in Python this is usually not what you want. 153 Do not supply the 'int', 'default', and 'maxwidth' arguments. 154 """ 155 156 # This code is a bit messy to make it fast for the 157 # common case while still doing adequate error checking. 158 istart = int(start) 159 if istart != start: 160 raise ValueError, "non-integer arg 1 for randrange()" 161 if stop is default: 162 if istart > 0: 163 if istart >= maxwidth: 164 return self._randbelow(istart) 165 return int(self.random() * istart) 166 raise ValueError, "empty range for randrange()" 167 168 # stop argument supplied. 169 istop = int(stop) 170 if istop != stop: 171 raise ValueError, "non-integer stop for randrange()" 172 width = istop - istart 173 if step == 1 and width > 0: 174 # Note that 175 # int(istart + self.random()*width) 176 # instead would be incorrect. For example, consider istart 177 # = -2 and istop = 0. Then the guts would be in 178 # -2.0 to 0.0 exclusive on both ends (ignoring that random() 179 # might return 0.0), and because int() truncates toward 0, the 180 # final result would be -1 or 0 (instead of -2 or -1). 181 # istart + int(self.random()*width) 182 # would also be incorrect, for a subtler reason: the RHS 183 # can return a long, and then randrange() would also return 184 # a long, but we're supposed to return an int (for backward 185 # compatibility). 186 187 if width >= maxwidth: 188 return int(istart + self._randbelow(width)) 189 return int(istart + int(self.random()*width)) 190 if step == 1: 191 raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width) 192 193 # Non-unit step argument supplied. 194 istep = int(step) 195 if istep != step: 196 raise ValueError, "non-integer step for randrange()" 197 if istep > 0: 198 n = (width + istep - 1) // istep 199 elif istep < 0: 200 n = (width + istep + 1) // istep 201 else: 202 raise ValueError, "zero step for randrange()" 203 204 if n <= 0: 205 raise ValueError, "empty range for randrange()" 206 207 if n >= maxwidth: 208 return istart + self._randbelow(n) 209 return istart + istep*int(self.random() * n)

210

211 - def randint(self, a, b):

212 """Return random integer in range [a, b], including both end points. 213 """ 214 215 return self.randrange(a, b+1)

216

217 - def _randbelow(self, n, _log=_log, int=int, _maxwidth=1L<<BPF, 218 _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType):

219 """Return a random int in the range [0,n) 220 221 Handles the case where n has more bits than returned 222 by a single call to the underlying generator. 223 """ 224 225 try: 226 getrandbits = self.getrandbits 227 except AttributeError: 228 pass 229 else: 230 # Only call self.getrandbits if the original random() builtin method 231 # has not been overridden or if a new getrandbits() was supplied. 232 # This assures that the two methods correspond. 233 if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method: 234 k = int(1.00001 + _log(n-1, 2.0)) # 2**k > n-1 > 2**(k-2) 235 r = getrandbits(k) 236 while r >= n: 237 r = getrandbits(k) 238 return r 239 if n >= _maxwidth: 240 _warn("Underlying random() generator does not supply \n" 241 "enough bits to choose from a population range this large") 242 return int(self.random() * n)

243 244 ## -------------------- sequence methods ------------------- 245

246 - def choice(self, seq):

247 """Choose a random element from a non-empty sequence.""" 248 return seq[int(self.random() * len(seq))] # raises IndexError if seq is empty

249

250 - def shuffle(self, x, random=None, int=int):

251 """x, random=random.random -> shuffle list x in place; return None. 252 253 Optional arg random is a 0-argument function returning a random 254 float in [0.0, 1.0); by default, the standard random.random. 255 """ 256 257 if random is None: 258 random = self.random 259 for i in reversed(xrange(1, len(x))): 260 # pick an element in x[:i+1] with which to exchange x[i] 261 j = int(random() * (i+1)) 262 x[i], x[j] = x[j], x[i]

263

264 - def sample(self, population, k):

265 """Chooses k unique random elements from a population sequence. 266 267 Returns a new list containing elements from the population while 268 leaving the original population unchanged. The resulting list is 269 in selection order so that all sub-slices will also be valid random 270 samples. This allows raffle winners (the sample) to be partitioned 271 into grand prize and second place winners (the subslices). 272 273 Members of the population need not be hashable or unique. If the 274 population contains repeats, then each occurrence is a possible 275 selection in the sample. 276 277 To choose a sample in a range of integers, use xrange as an argument. 278 This is especially fast and space efficient for sampling from a 279 large population: sample(xrange(10000000), 60) 280 """ 281 282 # XXX Although the documentation says `population` is "a sequence", 283 # XXX attempts are made to cater to any iterable with a __len__ 284 # XXX method. This has had mixed success. Examples from both 285 # XXX sides: sets work fine, and should become officially supported; 286 # XXX dicts are much harder, and have failed in various subtle 287 # XXX ways across attempts. Support for mapping types should probably 288 # XXX be dropped (and users should pass mapping.keys() or .values() 289 # XXX explicitly). 290 291 # Sampling without replacement entails tracking either potential 292 # selections (the pool) in a list or previous selections in a 293 # dictionary. 294 295 # When the number of selections is small compared to the 296 # population, then tracking selections is efficient, requiring 297 # only a small dictionary and an occasional reselection. For 298 # a larger number of selections, the pool tracking method is 299 # preferred since the list takes less space than the 300 # dictionary and it doesn't suffer from frequent reselections. 301 302 n = len(population) 303 if not 0 <= k <= n: 304 raise ValueError, "sample larger than population" 305 random = self.random 306 _int = int 307 result = [None] * k 308 if n < 6 * k or hasattr(population, "keys"): 309 # An n-length list is smaller than a k-length set, or this is a 310 # mapping type so the other algorithm wouldn't work. 311 pool = list(population) 312 for i in xrange(k): # invariant: non-selected at [0,n-i) 313 j = _int(random() * (n-i)) 314 result[i] = pool[j] 315 pool[j] = pool[n-i-1] # move non-selected item into vacancy 316 else: 317 try: 318 selected = {} 319 for i in xrange(k): 320 j = _int(random() * n) 321 while j in selected: 322 j = _int(random() * n) 323 result[i] = selected[j] = population[j] 324 except (TypeError, KeyError): # handle (at least) sets 325 if isinstance(population, list): 326 raise 327 return self.sample(tuple(population), k) 328 return result

329 330 ## -------------------- real-valued distributions ------------------- 331 332 ## -------------------- uniform distribution ------------------- 333

334 - def uniform(self, a, b):

335 """Get a random number in the range [a, b).""" 336 return a + (b-a) * self.random()

337 338 ## -------------------- normal distribution -------------------- 339

340 - def normalvariate(self, mu, sigma):

341 """Normal distribution. 342 343 mu is the mean, and sigma is the standard deviation. 344 345 """ 346 # mu = mean, sigma = standard deviation 347 348 # Uses Kinderman and Monahan method. Reference: Kinderman, 349 # A.J. and Monahan, J.F., "Computer generation of random 350 # variables using the ratio of uniform deviates", ACM Trans 351 # Math Software, 3, (1977), pp257-260. 352 353 random = self.random 354 while 1: 355 u1 = random() 356 u2 = 1.0 - random() 357 z = NV_MAGICCONST*(u1-0.5)/u2 358 zz = z*z/4.0 359 if zz <= -_log(u2): 360 break 361 return mu + z*sigma

362 363 ## -------------------- lognormal distribution -------------------- 364

365 - def lognormvariate(self, mu, sigma):

366 """Log normal distribution. 367 368 If you take the natural logarithm of this distribution, you'll get a 369 normal distribution with mean mu and standard deviation sigma. 370 mu can have any value, and sigma must be greater than zero. 371 372 """ 373 return _exp(self.normalvariate(mu, sigma))

374 375 ## -------------------- exponential distribution -------------------- 376

377 - def expovariate(self, lambd):

378 """Exponential distribution. 379 380 lambd is 1.0 divided by the desired mean. (The parameter would be 381 called "lambda", but that is a reserved word in Python.) Returned 382 values range from 0 to positive infinity. 383 384 """ 385 # lambd: rate lambd = 1/mean 386 # ('lambda' is a Python reserved word) 387 388 random = self.random 389 u = random() 390 while u <= 1e-7: 391 u = random() 392 return -_log(u)/lambd

393 394 ## -------------------- von Mises distribution -------------------- 395

396 - def vonmisesvariate(self, mu, kappa):

397 """Circular data distribution. 398 399 mu is the mean angle, expressed in radians between 0 and 2*pi, and 400 kappa is the concentration parameter, which must be greater than or 401 equal to zero. If kappa is equal to zero, this distribution reduces 402 to a uniform random angle over the range 0 to 2*pi. 403 404 """ 405 # mu: mean angle (in radians between 0 and 2*pi) 406 # kappa: concentration parameter kappa (>= 0) 407 # if kappa = 0 generate uniform random angle 408 409 # Based upon an algorithm published in: Fisher, N.I., 410 # "Statistical Analysis of Circular Data", Cambridge 411 # University Press, 1993. 412 413 # Thanks to Magnus Kessler for a correction to the 414 # implementation of step 4. 415 416 random = self.random 417 if kappa <= 1e-6: 418 return TWOPI * random() 419 420 a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa) 421 b = (a - _sqrt(2.0 * a))/(2.0 * kappa) 422 r = (1.0 + b * b)/(2.0 * b) 423 424 while 1: 425 u1 = random() 426 427 z = _cos(_pi * u1) 428 f = (1.0 + r * z)/(r + z) 429 c = kappa * (r - f) 430 431 u2 = random() 432 433 if u2 < c * (2.0 - c) or u2 <= c * _exp(1.0 - c): 434 break 435 436 u3 = random() 437 if u3 > 0.5: 438 theta = (mu % TWOPI) + _acos(f) 439 else: 440 theta = (mu % TWOPI) - _acos(f) 441 442 return theta

443 444 ## -------------------- gamma distribution -------------------- 445

446 - def gammavariate(self, alpha, beta):

447 """Gamma distribution. Not the gamma function! 448 449 Conditions on the parameters are alpha > 0 and beta > 0. 450 451 """ 452 453 # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 454 455 # Warning: a few older sources define the gamma distribution in terms 456 # of alpha > -1.0 457 if alpha <= 0.0 or beta <= 0.0: 458 raise ValueError, 'gammavariate: alpha and beta must be > 0.0' 459 460 random = self.random 461 if alpha > 1.0: 462 463 # Uses R.C.H. Cheng, "The generation of Gamma 464 # variables with non-integral shape parameters", 465 # Applied Statistics, (1977), 26, No. 1, p71-74 466 467 ainv = _sqrt(2.0 * alpha - 1.0) 468 bbb = alpha - LOG4 469 ccc = alpha + ainv 470 471 while 1: 472 u1 = random() 473 if not 1e-7 < u1 < .9999999: 474 continue 475 u2 = 1.0 - random() 476 v = _log(u1/(1.0-u1))/ainv 477 x = alpha*_exp(v) 478 z = u1*u1*u2 479 r = bbb+ccc*v-x 480 if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 481 return x * beta 482 483 elif alpha == 1.0: 484 # expovariate(1) 485 u = random() 486 while u <= 1e-7: 487 u = random() 488 return -_log(u) * beta 489 490 else: # alpha is between 0 and 1 (exclusive) 491 492 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 493 494 while 1: 495 u = random() 496 b = (_e + alpha)/_e 497 p = b*u 498 if p <= 1.0: 499 x = p ** (1.0/alpha) 500 else: 501 x = -_log((b-p)/alpha) 502 u1 = random() 503 if p > 1.0: 504 if u1 <= x ** (alpha - 1.0): 505 break 506 elif u1 <= _exp(-x): 507 break 508 return x * beta

509 510 ## -------------------- Gauss (faster alternative) -------------------- 511

512 - def gauss(self, mu, sigma):

513 """Gaussian distribution. 514 515 mu is the mean, and sigma is the standard deviation. This is 516 slightly faster than the normalvariate() function. 517 518 Not thread-safe without a lock around calls. 519 520 """ 521 522 # When x and y are two variables from [0, 1), uniformly 523 # distributed, then 524 # 525 # cos(2*pi*x)*sqrt(-2*log(1-y)) 526 # sin(2*pi*x)*sqrt(-2*log(1-y)) 527 # 528 # are two *independent* variables with normal distribution 529 # (mu = 0, sigma = 1). 530 # (Lambert Meertens) 531 # (corrected version; bug discovered by Mike Miller, fixed by LM) 532 533 # Multithreading note: When two threads call this function 534 # simultaneously, it is possible that they will receive the 535 # same return value. The window is very small though. To 536 # avoid this, you have to use a lock around all calls. (I 537 # didn't want to slow this down in the serial case by using a 538 # lock here.) 539 540 random = self.random 541 z = self.gauss_next 542 self.gauss_next = None 543 if z is None: 544 x2pi = random() * TWOPI 545 g2rad = _sqrt(-2.0 * _log(1.0 - random())) 546 z = _cos(x2pi) * g2rad 547 self.gauss_next = _sin(x2pi) * g2rad 548 549 return mu + z*sigma

550 551 ## -------------------- beta -------------------- 552 ## See 553 ## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470 554 ## for Ivan Frohne's insightful analysis of why the original implementation: 555 ## 556 ## def betavariate(self, alpha, beta): 557 ## # Discrete Event Simulation in C, pp 87-88. 558 ## 559 ## y = self.expovariate(alpha) 560 ## z = self.expovariate(1.0/beta) 561 ## return z/(y+z) 562 ## 563 ## was dead wrong, and how it probably got that way. 564

565 - def betavariate(self, alpha, beta):

566 """Beta distribution. 567 568 Conditions on the parameters are alpha > -1 and beta} > -1. 569 Returned values range between 0 and 1. 570 571 """ 572 573 # This version due to Janne Sinkkonen, and matches all the std 574 # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 575 y = self.gammavariate(alpha, 1.) 576 if y == 0: 577 return 0.0 578 else: 579 return y / (y + self.gammavariate(beta, 1.))

580 581 ## -------------------- Pareto -------------------- 582

583 - def paretovariate(self, alpha):

584 """Pareto distribution. alpha is the shape parameter.""" 585 # Jain, pg. 495 586 587 u = 1.0 - self.random() 588 return 1.0 / pow(u, 1.0/alpha)

589 590 ## -------------------- Weibull -------------------- 591

592 - def weibullvariate(self, alpha, beta):

593 """Weibull distribution. 594 595 alpha is the scale parameter and beta is the shape parameter. 596 597 """ 598 # Jain, pg. 499; bug fix courtesy Bill Arms 599 600 u = 1.0 - self.random() 601 return alpha * pow(-_log(u), 1.0/beta)

602 603 ## -------------------- Wichmann-Hill ------------------- 604

605 -class WichmannHill(Random):

606 607 VERSION = 1 # used by getstate/setstate 608

609 - def seed(self, a=None):

610 """Initialize internal state from hashable object. 611 612 None or no argument seeds from current time or from an operating 613 system specific randomness source if available. 614 615 If a is not None or an int or long, hash(a) is used instead. 616 617 If a is an int or long, a is used directly. Distinct values between 618 0 and 27814431486575L inclusive are guaranteed to yield distinct 619 internal states (this guarantee is specific to the default 620 Wichmann-Hill generator). 621 """ 622 623 if a is None: 624 try: 625 a = long(_hexlify(_urandom(16)), 16) 626 except NotImplementedError: 627 import time 628 a = long(time.time() * 256) # use fractional seconds 629 630 if not isinstance(a, (int, long)): 631 a = hash(a) 632 633 a, x = divmod(a, 30268) 634 a, y = divmod(a, 30306) 635 a, z = divmod(a, 30322) 636 self._seed = int(x)+1, int(y)+1, int(z)+1 637 638 self.gauss_next = None

639

640 - def random(self):

641 """Get the next random number in the range [0.0, 1.0).""" 642 643 # Wichman-Hill random number generator. 644 # 645 # Wichmann, B. A. & Hill, I. D. (1982) 646 # Algorithm AS 183: 647 # An efficient and portable pseudo-random number generator 648 # Applied Statistics 31 (1982) 188-190 649 # 650 # see also: 651 # Correction to Algorithm AS 183 652 # Applied Statistics 33 (1984) 123 653 # 654 # McLeod, A. I. (1985) 655 # A remark on Algorithm AS 183 656 # Applied Statistics 34 (1985),198-200 657 658 # This part is thread-unsafe: 659 # BEGIN CRITICAL SECTION 660 x, y, z = self._seed 661 x = (171 * x) % 30269 662 y = (172 * y) % 30307 663 z = (170 * z) % 30323 664 self._seed = x, y, z 665 # END CRITICAL SECTION 666 667 # Note: on a platform using IEEE-754 double arithmetic, this can 668 # never return 0.0 (asserted by Tim; proof too long for a comment). 669 return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0

670

671 - def getstate(self):

672 """Return internal state; can be passed to setstate() later.""" 673 return self.VERSION, self._seed, self.gauss_next

674

675 - def setstate(self, state):

676 """Restore internal state from object returned by getstate().""" 677 version = state[0] 678 if version == 1: 679 version, self._seed, self.gauss_next = state 680 else: 681 raise ValueError("state with version %s passed to " 682 "Random.setstate() of version %s" % 683 (version, self.VERSION))

684

685 - def jumpahead(self, n):

686 """Act as if n calls to random() were made, but quickly. 687 688 n is an int, greater than or equal to 0. 689 690 Example use: If you have 2 threads and know that each will 691 consume no more than a million random numbers, create two Random 692 objects r1 and r2, then do 693 r2.setstate(r1.getstate()) 694 r2.jumpahead(1000000) 695 Then r1 and r2 will use guaranteed-disjoint segments of the full 696 period. 697 """ 698 699 if not n >= 0: 700 raise ValueError("n must be >= 0") 701 x, y, z = self._seed 702 x = int(x * pow(171, n, 30269)) % 30269 703 y = int(y * pow(172, n, 30307)) % 30307 704 z = int(z * pow(170, n, 30323)) % 30323 705 self._seed = x, y, z

706

707 - def __whseed(self, x=0, y=0, z=0):

708 """Set the Wichmann-Hill seed from (x, y, z). 709 710 These must be integers in the range [0, 256). 711 """ 712 713 if not type(x) == type(y) == type(z) == int: 714 raise TypeError('seeds must be integers') 715 if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): 716 raise ValueError('seeds must be in range(0, 256)') 717 if 0 == x == y == z: 718 # Initialize from current time 719 import time 720 t = long(time.time() * 256) 721 t = int((t&0xffffff) ^ (t>>24)) 722 t, x = divmod(t, 256) 723 t, y = divmod(t, 256) 724 t, z = divmod(t, 256) 725 # Zero is a poor seed, so substitute 1 726 self._seed = (x or 1, y or 1, z or 1) 727 728 self.gauss_next = None

729

730 - def whseed(self, a=None):

731 """Seed from hashable object's hash code. 732 733 None or no argument seeds from current time. It is not guaranteed 734 that objects with distinct hash codes lead to distinct internal 735 states. 736 737 This is obsolete, provided for compatibility with the seed routine 738 used prior to Python 2.1. Use the .seed() method instead. 739 """ 740 741 if a is None: 742 self.__whseed() 743 return 744 a = hash(a) 745 a, x = divmod(a, 256) 746 a, y = divmod(a, 256) 747 a, z = divmod(a, 256) 748 x = (x + a) % 256 or 1 749 y = (y + a) % 256 or 1 750 z = (z + a) % 256 or 1 751 self.__whseed(x, y, z)

752 753 ## --------------- Operating System Random Source ------------------ 754

755 -class SystemRandom(Random):

756 """Alternate random number generator using sources provided 757 by the operating system (such as /dev/urandom on Unix or 758 CryptGenRandom on Windows). 759 760 Not available on all systems (see os.urandom() for details). 761 """ 762

763 - def random(self):

764 """Get the next random number in the range [0.0, 1.0).""" 765 return (long(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF

766

767 - def getrandbits(self, k):

768 """getrandbits(k) -> x. Generates a long int with k random bits.""" 769 if k <= 0: 770 raise ValueError('number of bits must be greater than zero') 771 if k != int(k): 772 raise TypeError('number of bits should be an integer') 773 bytes = (k + 7) // 8 # bits / 8 and rounded up 774 x = long(_hexlify(_urandom(bytes)), 16) 775 return x >> (bytes * 8 - k) # trim excess bits

776

777 - def _stub(self, *args, **kwds):

778 "Stub method. Not used for a system random number generator." 779 return None

780 seed = jumpahead = _stub 781

782 - def _notimplemented(self, *args, **kwds):

783 "Method should not be called for a system random number generator." 784 raise NotImplementedError('System entropy source does not have state.')

785 getstate = setstate = _notimplemented

786 787 ## -------------------- test program -------------------- 788

789 -def _test_generator(n, func, args):

790 import time 791 print n, 'times', func.__name__ 792 total = 0.0 793 sqsum = 0.0 794 smallest = 1e10 795 largest = -1e10 796 t0 = time.time() 797 for i in range(n): 798 x = func(*args) 799 total += x 800 sqsum = sqsum + x*x 801 smallest = min(x, smallest) 802 largest = max(x, largest) 803 t1 = time.time() 804 print round(t1-t0, 3), 'sec,', 805 avg = total/n 806 stddev = _sqrt(sqsum/n - avg*avg) 807 print 'avg %g, stddev %g, min %g, max %g' % \ 808 (avg, stddev, smallest, largest)

809 810

811 -def _test(N=2000):

812 _test_generator(N, random, ()) 813 _test_generator(N, normalvariate, (0.0, 1.0)) 814 _test_generator(N, lognormvariate, (0.0, 1.0)) 815 _test_generator(N, vonmisesvariate, (0.0, 1.0)) 816 _test_generator(N, gammavariate, (0.01, 1.0)) 817 _test_generator(N, gammavariate, (0.1, 1.0)) 818 _test_generator(N, gammavariate, (0.1, 2.0)) 819 _test_generator(N, gammavariate, (0.5, 1.0)) 820 _test_generator(N, gammavariate, (0.9, 1.0)) 821 _test_generator(N, gammavariate, (1.0, 1.0)) 822 _test_generator(N, gammavariate, (2.0, 1.0)) 823 _test_generator(N, gammavariate, (20.0, 1.0)) 824 _test_generator(N, gammavariate, (200.0, 1.0)) 825 _test_generator(N, gauss, (0.0, 1.0)) 826 _test_generator(N, betavariate, (3.0, 3.0))

827 828 # Create one instance, seeded from current time, and export its methods 829 # as module-level functions. The functions share state across all uses 830 #(both in the user's code and in the Python libraries), but that's fine 831 # for most programs and is easier for the casual user than making them 832 # instantiate their own Random() instance. 833 834 _inst = Random() 835 seed = _inst.seed 836 random = _inst.random 837 uniform = _inst.uniform 838 randint = _inst.randint 839 choice = _inst.choice 840 randrange = _inst.randrange 841 sample = _inst.sample 842 shuffle = _inst.shuffle 843 normalvariate = _inst.normalvariate 844 lognormvariate = _inst.lognormvariate 845 expovariate = _inst.expovariate 846 vonmisesvariate = _inst.vonmisesvariate 847 gammavariate = _inst.gammavariate 848 gauss = _inst.gauss 849 betavariate = _inst.betavariate 850 paretovariate = _inst.paretovariate 851 weibullvariate = _inst.weibullvariate 852 getstate = _inst.getstate 853 setstate = _inst.setstate 854 jumpahead = _inst.jumpahead 855 getrandbits = _inst.getrandbits 856 857 if __name__ == '__main__': 858 _test() 859

Source Code for Module random