random.py revision c0b4034b8165ce958a23f2c865b51ae0f52040f5
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 22 distributions on the circle (angles 0 to 2pi) 23 --------------------------------------------- 24 circular uniform 25 von Mises 26 27Translated from anonymously contributed C/C++ source. 28 29Multi-threading note: the random number generator used here is not thread- 30safe; it is possible that two calls return the same random value. However, 31you can instantiate a different instance of Random() in each thread to get 32generators that don't share state, then use .setstate() and .jumpahead() to 33move the generators to disjoint segments of the full period. For example, 34 35def create_generators(num, delta, firstseed=None): 36 ""\"Return list of num distinct generators. 37 Each generator has its own unique segment of delta elements from 38 Random.random()'s full period. 39 Seed the first generator with optional arg firstseed (default is 40 None, to seed from current time). 41 ""\" 42 43 from random import Random 44 g = Random(firstseed) 45 result = [g] 46 for i in range(num - 1): 47 laststate = g.getstate() 48 g = Random() 49 g.setstate(laststate) 50 g.jumpahead(delta) 51 result.append(g) 52 return result 53 54gens = create_generators(10, 1000000) 55 56That creates 10 distinct generators, which can be passed out to 10 distinct 57threads. The generators don't share state so can be called safely in 58parallel. So long as no thread calls its g.random() more than a million 59times (the second argument to create_generators), the sequences seen by 60each thread will not overlap. 61 62The period of the underlying Wichmann-Hill generator is 6,953,607,871,644, 63and that limits how far this technique can be pushed. 64 65Just for fun, note that since we know the period, .jumpahead() can also be 66used to "move backward in time": 67 68>>> g = Random(42) # arbitrary 69>>> g.random() 700.25420336316883324 71>>> g.jumpahead(6953607871644L - 1) # move *back* one 72>>> g.random() 730.25420336316883324 74""" 75# XXX The docstring sucks. 76 77from math import log as _log, exp as _exp, pi as _pi, e as _e 78from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin 79from math import floor as _floor 80 81__all__ = ["Random","seed","random","uniform","randint","choice","sample", 82 "randrange","shuffle","normalvariate","lognormvariate", 83 "cunifvariate","expovariate","vonmisesvariate","gammavariate", 84 "stdgamma","gauss","betavariate","paretovariate","weibullvariate", 85 "getstate","setstate","jumpahead","whseed"] 86 87def _verify(name, computed, expected): 88 if abs(computed - expected) > 1e-7: 89 raise ValueError( 90 "computed value for %s deviates too much " 91 "(computed %g, expected %g)" % (name, computed, expected)) 92 93NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) 94_verify('NV_MAGICCONST', NV_MAGICCONST, 1.71552776992141) 95 96TWOPI = 2.0*_pi 97_verify('TWOPI', TWOPI, 6.28318530718) 98 99LOG4 = _log(4.0) 100_verify('LOG4', LOG4, 1.38629436111989) 101 102SG_MAGICCONST = 1.0 + _log(4.5) 103_verify('SG_MAGICCONST', SG_MAGICCONST, 2.50407739677627) 104 105del _verify 106 107# Translated by Guido van Rossum from C source provided by 108# Adrian Baddeley. 109 110class Random: 111 """Random number generator base class used by bound module functions. 112 113 Used to instantiate instances of Random to get generators that don't 114 share state. Especially useful for multi-threaded programs, creating 115 a different instance of Random for each thread, and using the jumpahead() 116 method to ensure that the generated sequences seen by each thread don't 117 overlap. 118 119 Class Random can also be subclassed if you want to use a different basic 120 generator of your own devising: in that case, override the following 121 methods: random(), seed(), getstate(), setstate() and jumpahead(). 122 123 """ 124 125 VERSION = 1 # used by getstate/setstate 126 127 def __init__(self, x=None): 128 """Initialize an instance. 129 130 Optional argument x controls seeding, as for Random.seed(). 131 """ 132 133 self.seed(x) 134 135## -------------------- core generator ------------------- 136 137 # Specific to Wichmann-Hill generator. Subclasses wishing to use a 138 # different core generator should override the seed(), random(), 139 # getstate(), setstate() and jumpahead() methods. 140 141 def seed(self, a=None): 142 """Initialize internal state from hashable object. 143 144 None or no argument seeds from current time. 145 146 If a is not None or an int or long, hash(a) is used instead. 147 148 If a is an int or long, a is used directly. Distinct values between 149 0 and 27814431486575L inclusive are guaranteed to yield distinct 150 internal states (this guarantee is specific to the default 151 Wichmann-Hill generator). 152 """ 153 154 if a is None: 155 # Initialize from current time 156 import time 157 a = long(time.time() * 256) 158 159 if type(a) not in (type(3), type(3L)): 160 a = hash(a) 161 162 a, x = divmod(a, 30268) 163 a, y = divmod(a, 30306) 164 a, z = divmod(a, 30322) 165 self._seed = int(x)+1, int(y)+1, int(z)+1 166 167 self.gauss_next = None 168 169 def random(self): 170 """Get the next random number in the range [0.0, 1.0).""" 171 172 # Wichman-Hill random number generator. 173 # 174 # Wichmann, B. A. & Hill, I. D. (1982) 175 # Algorithm AS 183: 176 # An efficient and portable pseudo-random number generator 177 # Applied Statistics 31 (1982) 188-190 178 # 179 # see also: 180 # Correction to Algorithm AS 183 181 # Applied Statistics 33 (1984) 123 182 # 183 # McLeod, A. I. (1985) 184 # A remark on Algorithm AS 183 185 # Applied Statistics 34 (1985),198-200 186 187 # This part is thread-unsafe: 188 # BEGIN CRITICAL SECTION 189 x, y, z = self._seed 190 x = (171 * x) % 30269 191 y = (172 * y) % 30307 192 z = (170 * z) % 30323 193 self._seed = x, y, z 194 # END CRITICAL SECTION 195 196 # Note: on a platform using IEEE-754 double arithmetic, this can 197 # never return 0.0 (asserted by Tim; proof too long for a comment). 198 return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0 199 200 def getstate(self): 201 """Return internal state; can be passed to setstate() later.""" 202 return self.VERSION, self._seed, self.gauss_next 203 204 def setstate(self, state): 205 """Restore internal state from object returned by getstate().""" 206 version = state[0] 207 if version == 1: 208 version, self._seed, self.gauss_next = state 209 else: 210 raise ValueError("state with version %s passed to " 211 "Random.setstate() of version %s" % 212 (version, self.VERSION)) 213 214 def jumpahead(self, n): 215 """Act as if n calls to random() were made, but quickly. 216 217 n is an int, greater than or equal to 0. 218 219 Example use: If you have 2 threads and know that each will 220 consume no more than a million random numbers, create two Random 221 objects r1 and r2, then do 222 r2.setstate(r1.getstate()) 223 r2.jumpahead(1000000) 224 Then r1 and r2 will use guaranteed-disjoint segments of the full 225 period. 226 """ 227 228 if not n >= 0: 229 raise ValueError("n must be >= 0") 230 x, y, z = self._seed 231 x = int(x * pow(171, n, 30269)) % 30269 232 y = int(y * pow(172, n, 30307)) % 30307 233 z = int(z * pow(170, n, 30323)) % 30323 234 self._seed = x, y, z 235 236 def __whseed(self, x=0, y=0, z=0): 237 """Set the Wichmann-Hill seed from (x, y, z). 238 239 These must be integers in the range [0, 256). 240 """ 241 242 if not type(x) == type(y) == type(z) == type(0): 243 raise TypeError('seeds must be integers') 244 if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): 245 raise ValueError('seeds must be in range(0, 256)') 246 if 0 == x == y == z: 247 # Initialize from current time 248 import time 249 t = long(time.time() * 256) 250 t = int((t&0xffffff) ^ (t>>24)) 251 t, x = divmod(t, 256) 252 t, y = divmod(t, 256) 253 t, z = divmod(t, 256) 254 # Zero is a poor seed, so substitute 1 255 self._seed = (x or 1, y or 1, z or 1) 256 257 self.gauss_next = None 258 259 def whseed(self, a=None): 260 """Seed from hashable object's hash code. 261 262 None or no argument seeds from current time. It is not guaranteed 263 that objects with distinct hash codes lead to distinct internal 264 states. 265 266 This is obsolete, provided for compatibility with the seed routine 267 used prior to Python 2.1. Use the .seed() method instead. 268 """ 269 270 if a is None: 271 self.__whseed() 272 return 273 a = hash(a) 274 a, x = divmod(a, 256) 275 a, y = divmod(a, 256) 276 a, z = divmod(a, 256) 277 x = (x + a) % 256 or 1 278 y = (y + a) % 256 or 1 279 z = (z + a) % 256 or 1 280 self.__whseed(x, y, z) 281 282## ---- Methods below this point do not need to be overridden when 283## ---- subclassing for the purpose of using a different core generator. 284 285## -------------------- pickle support ------------------- 286 287 def __getstate__(self): # for pickle 288 return self.getstate() 289 290 def __setstate__(self, state): # for pickle 291 self.setstate(state) 292 293## -------------------- integer methods ------------------- 294 295 def randrange(self, start, stop=None, step=1, int=int, default=None): 296 """Choose a random item from range(start, stop[, step]). 297 298 This fixes the problem with randint() which includes the 299 endpoint; in Python this is usually not what you want. 300 Do not supply the 'int' and 'default' arguments. 301 """ 302 303 # This code is a bit messy to make it fast for the 304 # common case while still doing adequate error checking. 305 istart = int(start) 306 if istart != start: 307 raise ValueError, "non-integer arg 1 for randrange()" 308 if stop is default: 309 if istart > 0: 310 return int(self.random() * istart) 311 raise ValueError, "empty range for randrange()" 312 313 # stop argument supplied. 314 istop = int(stop) 315 if istop != stop: 316 raise ValueError, "non-integer stop for randrange()" 317 if step == 1 and istart < istop: 318 try: 319 return istart + int(self.random()*(istop - istart)) 320 except OverflowError: 321 # This can happen if istop-istart > sys.maxint + 1, and 322 # multiplying by random() doesn't reduce it to something 323 # <= sys.maxint. We know that the overall result fits 324 # in an int, and can still do it correctly via math.floor(). 325 # But that adds another function call, so for speed we 326 # avoided that whenever possible. 327 return int(istart + _floor(self.random()*(istop - istart))) 328 if step == 1: 329 raise ValueError, "empty range for randrange()" 330 331 # Non-unit step argument supplied. 332 istep = int(step) 333 if istep != step: 334 raise ValueError, "non-integer step for randrange()" 335 if istep > 0: 336 n = (istop - istart + istep - 1) / istep 337 elif istep < 0: 338 n = (istop - istart + istep + 1) / istep 339 else: 340 raise ValueError, "zero step for randrange()" 341 342 if n <= 0: 343 raise ValueError, "empty range for randrange()" 344 return istart + istep*int(self.random() * n) 345 346 def randint(self, a, b): 347 """Return random integer in range [a, b], including both end points. 348 """ 349 350 return self.randrange(a, b+1) 351 352## -------------------- sequence methods ------------------- 353 354 def choice(self, seq): 355 """Choose a random element from a non-empty sequence.""" 356 return seq[int(self.random() * len(seq))] 357 358 def shuffle(self, x, random=None, int=int): 359 """x, random=random.random -> shuffle list x in place; return None. 360 361 Optional arg random is a 0-argument function returning a random 362 float in [0.0, 1.0); by default, the standard random.random. 363 364 Note that for even rather small len(x), the total number of 365 permutations of x is larger than the period of most random number 366 generators; this implies that "most" permutations of a long 367 sequence can never be generated. 368 """ 369 370 if random is None: 371 random = self.random 372 for i in xrange(len(x)-1, 0, -1): 373 # pick an element in x[:i+1] with which to exchange x[i] 374 j = int(random() * (i+1)) 375 x[i], x[j] = x[j], x[i] 376 377 def sample(self, population, k, random=None, int=int): 378 """Chooses k unique random elements from a population sequence. 379 380 Returns a new list containing elements from the population while 381 leaving the original population unchanged. The resulting list is 382 in selection order so that all sub-slices will also be valid random 383 samples. This allows raffle winners (the sample) to be partitioned 384 into grand prize and second place winners (the subslices). 385 386 Members of the population need not be hashable or unique. If the 387 population contains repeats, then each occurrence is a possible 388 selection in the sample. 389 390 To choose a sample in a range of integers, use xrange as an argument. 391 This is especially fast and space efficient for sampling from a 392 large population: sample(xrange(10000000), 60) 393 394 Optional arg random is a 0-argument function returning a random 395 float in [0.0, 1.0); by default, the standard random.random. 396 """ 397 398 # Sampling without replacement entails tracking either potential 399 # selections (the pool) or previous selections. 400 401 # Pools are stored in lists which provide __getitem__ for selection 402 # and provide a way to remove selections. But each list.remove() 403 # rebuilds the entire list, so it is better to rearrange the list, 404 # placing non-selected elements at the head of the list. Tracking 405 # the selection pool is only space efficient with small populations. 406 407 # Previous selections are stored in dictionaries which provide 408 # __contains__ for detecting repeat selections. Discarding repeats 409 # is efficient unless most of the population has already been chosen. 410 # So, tracking selections is useful when sample sizes are much 411 # smaller than the total population. 412 413 n = len(population) 414 if not 0 <= k <= n: 415 raise ValueError, "sample larger than population" 416 if random is None: 417 random = self.random 418 result = [None] * k 419 if n < 6 * k: # if n len list takes less space than a k len dict 420 pool = list(population) # track potential selections 421 for i in xrange(k): 422 j = int(random() * (n-i)) # non-selected at [0,n-i) 423 result[i] = pool[j] # save selected element 424 pool[j] = pool[n-i-1] # non-selected to head of list 425 else: 426 selected = {} # track previous selections 427 for i in xrange(k): 428 j = int(random() * n) 429 while j in selected: # discard and replace repeats 430 j = int(random() * n) 431 result[i] = selected[j] = population[j] 432 return result # return selections in the order they were picked 433 434## -------------------- real-valued distributions ------------------- 435 436## -------------------- uniform distribution ------------------- 437 438 def uniform(self, a, b): 439 """Get a random number in the range [a, b).""" 440 return a + (b-a) * self.random() 441 442## -------------------- normal distribution -------------------- 443 444 def normalvariate(self, mu, sigma): 445 """Normal distribution. 446 447 mu is the mean, and sigma is the standard deviation. 448 449 """ 450 # mu = mean, sigma = standard deviation 451 452 # Uses Kinderman and Monahan method. Reference: Kinderman, 453 # A.J. and Monahan, J.F., "Computer generation of random 454 # variables using the ratio of uniform deviates", ACM Trans 455 # Math Software, 3, (1977), pp257-260. 456 457 random = self.random 458 while 1: 459 u1 = random() 460 u2 = random() 461 z = NV_MAGICCONST*(u1-0.5)/u2 462 zz = z*z/4.0 463 if zz <= -_log(u2): 464 break 465 return mu + z*sigma 466 467## -------------------- lognormal distribution -------------------- 468 469 def lognormvariate(self, mu, sigma): 470 """Log normal distribution. 471 472 If you take the natural logarithm of this distribution, you'll get a 473 normal distribution with mean mu and standard deviation sigma. 474 mu can have any value, and sigma must be greater than zero. 475 476 """ 477 return _exp(self.normalvariate(mu, sigma)) 478 479## -------------------- circular uniform -------------------- 480 481 def cunifvariate(self, mean, arc): 482 """Circular uniform distribution. 483 484 mean is the mean angle, and arc is the range of the distribution, 485 centered around the mean angle. Both values must be expressed in 486 radians. Returned values range between mean - arc/2 and 487 mean + arc/2 and are normalized to between 0 and pi. 488 489 Deprecated in version 2.3. Use: 490 (mean + arc * (Random.random() - 0.5)) % Math.pi 491 492 """ 493 # mean: mean angle (in radians between 0 and pi) 494 # arc: range of distribution (in radians between 0 and pi) 495 import warnings 496 warnings.warn("The cunifvariate function is deprecated; Use (mean " 497 "+ arc * (Random.random() - 0.5)) % Math.pi instead", 498 DeprecationWarning) 499 500 return (mean + arc * (self.random() - 0.5)) % _pi 501 502## -------------------- exponential distribution -------------------- 503 504 def expovariate(self, lambd): 505 """Exponential distribution. 506 507 lambd is 1.0 divided by the desired mean. (The parameter would be 508 called "lambda", but that is a reserved word in Python.) Returned 509 values range from 0 to positive infinity. 510 511 """ 512 # lambd: rate lambd = 1/mean 513 # ('lambda' is a Python reserved word) 514 515 random = self.random 516 u = random() 517 while u <= 1e-7: 518 u = random() 519 return -_log(u)/lambd 520 521## -------------------- von Mises distribution -------------------- 522 523 def vonmisesvariate(self, mu, kappa): 524 """Circular data distribution. 525 526 mu is the mean angle, expressed in radians between 0 and 2*pi, and 527 kappa is the concentration parameter, which must be greater than or 528 equal to zero. If kappa is equal to zero, this distribution reduces 529 to a uniform random angle over the range 0 to 2*pi. 530 531 """ 532 # mu: mean angle (in radians between 0 and 2*pi) 533 # kappa: concentration parameter kappa (>= 0) 534 # if kappa = 0 generate uniform random angle 535 536 # Based upon an algorithm published in: Fisher, N.I., 537 # "Statistical Analysis of Circular Data", Cambridge 538 # University Press, 1993. 539 540 # Thanks to Magnus Kessler for a correction to the 541 # implementation of step 4. 542 543 random = self.random 544 if kappa <= 1e-6: 545 return TWOPI * random() 546 547 a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa) 548 b = (a - _sqrt(2.0 * a))/(2.0 * kappa) 549 r = (1.0 + b * b)/(2.0 * b) 550 551 while 1: 552 u1 = random() 553 554 z = _cos(_pi * u1) 555 f = (1.0 + r * z)/(r + z) 556 c = kappa * (r - f) 557 558 u2 = random() 559 560 if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)): 561 break 562 563 u3 = random() 564 if u3 > 0.5: 565 theta = (mu % TWOPI) + _acos(f) 566 else: 567 theta = (mu % TWOPI) - _acos(f) 568 569 return theta 570 571## -------------------- gamma distribution -------------------- 572 573 def gammavariate(self, alpha, beta): 574 """Gamma distribution. Not the gamma function! 575 576 Conditions on the parameters are alpha > 0 and beta > 0. 577 578 """ 579 580 # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 581 582 # Warning: a few older sources define the gamma distribution in terms 583 # of alpha > -1.0 584 if alpha <= 0.0 or beta <= 0.0: 585 raise ValueError, 'gammavariate: alpha and beta must be > 0.0' 586 587 random = self.random 588 if alpha > 1.0: 589 590 # Uses R.C.H. Cheng, "The generation of Gamma 591 # variables with non-integral shape parameters", 592 # Applied Statistics, (1977), 26, No. 1, p71-74 593 594 ainv = _sqrt(2.0 * alpha - 1.0) 595 bbb = alpha - LOG4 596 ccc = alpha + ainv 597 598 while 1: 599 u1 = random() 600 u2 = random() 601 v = _log(u1/(1.0-u1))/ainv 602 x = alpha*_exp(v) 603 z = u1*u1*u2 604 r = bbb+ccc*v-x 605 if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 606 return x * beta 607 608 elif alpha == 1.0: 609 # expovariate(1) 610 u = random() 611 while u <= 1e-7: 612 u = random() 613 return -_log(u) * beta 614 615 else: # alpha is between 0 and 1 (exclusive) 616 617 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 618 619 while 1: 620 u = random() 621 b = (_e + alpha)/_e 622 p = b*u 623 if p <= 1.0: 624 x = pow(p, 1.0/alpha) 625 else: 626 # p > 1 627 x = -_log((b-p)/alpha) 628 u1 = random() 629 if not (((p <= 1.0) and (u1 > _exp(-x))) or 630 ((p > 1) and (u1 > pow(x, alpha - 1.0)))): 631 break 632 return x * beta 633 634 635 def stdgamma(self, alpha, ainv, bbb, ccc): 636 # This method was (and shall remain) undocumented. 637 # This method is deprecated 638 # for the following reasons: 639 # 1. Returns same as .gammavariate(alpha, 1.0) 640 # 2. Requires caller to provide 3 extra arguments 641 # that are functions of alpha anyway 642 # 3. Can't be used for alpha < 0.5 643 644 # ainv = sqrt(2 * alpha - 1) 645 # bbb = alpha - log(4) 646 # ccc = alpha + ainv 647 import warnings 648 warnings.warn("The stdgamma function is deprecated; " 649 "use gammavariate() instead", 650 DeprecationWarning) 651 return self.gammavariate(alpha, 1.0) 652 653 654 655## -------------------- Gauss (faster alternative) -------------------- 656 657 def gauss(self, mu, sigma): 658 """Gaussian distribution. 659 660 mu is the mean, and sigma is the standard deviation. This is 661 slightly faster than the normalvariate() function. 662 663 Not thread-safe without a lock around calls. 664 665 """ 666 667 # When x and y are two variables from [0, 1), uniformly 668 # distributed, then 669 # 670 # cos(2*pi*x)*sqrt(-2*log(1-y)) 671 # sin(2*pi*x)*sqrt(-2*log(1-y)) 672 # 673 # are two *independent* variables with normal distribution 674 # (mu = 0, sigma = 1). 675 # (Lambert Meertens) 676 # (corrected version; bug discovered by Mike Miller, fixed by LM) 677 678 # Multithreading note: When two threads call this function 679 # simultaneously, it is possible that they will receive the 680 # same return value. The window is very small though. To 681 # avoid this, you have to use a lock around all calls. (I 682 # didn't want to slow this down in the serial case by using a 683 # lock here.) 684 685 random = self.random 686 z = self.gauss_next 687 self.gauss_next = None 688 if z is None: 689 x2pi = random() * TWOPI 690 g2rad = _sqrt(-2.0 * _log(1.0 - random())) 691 z = _cos(x2pi) * g2rad 692 self.gauss_next = _sin(x2pi) * g2rad 693 694 return mu + z*sigma 695 696## -------------------- beta -------------------- 697## See 698## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470 699## for Ivan Frohne's insightful analysis of why the original implementation: 700## 701## def betavariate(self, alpha, beta): 702## # Discrete Event Simulation in C, pp 87-88. 703## 704## y = self.expovariate(alpha) 705## z = self.expovariate(1.0/beta) 706## return z/(y+z) 707## 708## was dead wrong, and how it probably got that way. 709 710 def betavariate(self, alpha, beta): 711 """Beta distribution. 712 713 Conditions on the parameters are alpha > -1 and beta} > -1. 714 Returned values range between 0 and 1. 715 716 """ 717 718 # This version due to Janne Sinkkonen, and matches all the std 719 # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 720 y = self.gammavariate(alpha, 1.) 721 if y == 0: 722 return 0.0 723 else: 724 return y / (y + self.gammavariate(beta, 1.)) 725 726## -------------------- Pareto -------------------- 727 728 def paretovariate(self, alpha): 729 """Pareto distribution. alpha is the shape parameter.""" 730 # Jain, pg. 495 731 732 u = self.random() 733 return 1.0 / pow(u, 1.0/alpha) 734 735## -------------------- Weibull -------------------- 736 737 def weibullvariate(self, alpha, beta): 738 """Weibull distribution. 739 740 alpha is the scale parameter and beta is the shape parameter. 741 742 """ 743 # Jain, pg. 499; bug fix courtesy Bill Arms 744 745 u = self.random() 746 return alpha * pow(-_log(u), 1.0/beta) 747 748## -------------------- test program -------------------- 749 750def _test_generator(n, funccall): 751 import time 752 print n, 'times', funccall 753 code = compile(funccall, funccall, 'eval') 754 sum = 0.0 755 sqsum = 0.0 756 smallest = 1e10 757 largest = -1e10 758 t0 = time.time() 759 for i in range(n): 760 x = eval(code) 761 sum = sum + x 762 sqsum = sqsum + x*x 763 smallest = min(x, smallest) 764 largest = max(x, largest) 765 t1 = time.time() 766 print round(t1-t0, 3), 'sec,', 767 avg = sum/n 768 stddev = _sqrt(sqsum/n - avg*avg) 769 print 'avg %g, stddev %g, min %g, max %g' % \ 770 (avg, stddev, smallest, largest) 771 772def _test_sample(n): 773 # For the entire allowable range of 0 <= k <= n, validate that 774 # the sample is of the correct length and contains only unique items 775 population = xrange(n) 776 for k in xrange(n+1): 777 s = sample(population, k) 778 assert len(dict([(elem,True) for elem in s])) == len(s) == k 779 assert None not in s 780 781def _sample_generator(n, k): 782 # Return a fixed element from the sample. Validates random ordering. 783 return sample(xrange(n), k)[k//2] 784 785def _test(N=2000): 786 print 'TWOPI =', TWOPI 787 print 'LOG4 =', LOG4 788 print 'NV_MAGICCONST =', NV_MAGICCONST 789 print 'SG_MAGICCONST =', SG_MAGICCONST 790 _test_generator(N, 'random()') 791 _test_generator(N, 'normalvariate(0.0, 1.0)') 792 _test_generator(N, 'lognormvariate(0.0, 1.0)') 793 _test_generator(N, 'cunifvariate(0.0, 1.0)') 794 _test_generator(N, 'expovariate(1.0)') 795 _test_generator(N, 'vonmisesvariate(0.0, 1.0)') 796 _test_generator(N, 'gammavariate(0.01, 1.0)') 797 _test_generator(N, 'gammavariate(0.1, 1.0)') 798 _test_generator(N, 'gammavariate(0.1, 2.0)') 799 _test_generator(N, 'gammavariate(0.5, 1.0)') 800 _test_generator(N, 'gammavariate(0.9, 1.0)') 801 _test_generator(N, 'gammavariate(1.0, 1.0)') 802 _test_generator(N, 'gammavariate(2.0, 1.0)') 803 _test_generator(N, 'gammavariate(20.0, 1.0)') 804 _test_generator(N, 'gammavariate(200.0, 1.0)') 805 _test_generator(N, 'gauss(0.0, 1.0)') 806 _test_generator(N, 'betavariate(3.0, 3.0)') 807 _test_generator(N, 'paretovariate(1.0)') 808 _test_generator(N, 'weibullvariate(1.0, 1.0)') 809 _test_generator(N, '_sample_generator(50, 5)') # expected s.d.: 14.4 810 _test_generator(N, '_sample_generator(50, 45)') # expected s.d.: 14.4 811 _test_sample(500) 812 813 # Test jumpahead. 814 s = getstate() 815 jumpahead(N) 816 r1 = random() 817 # now do it the slow way 818 setstate(s) 819 for i in range(N): 820 random() 821 r2 = random() 822 if r1 != r2: 823 raise ValueError("jumpahead test failed " + `(N, r1, r2)`) 824 825# Create one instance, seeded from current time, and export its methods 826# as module-level functions. The functions are not threadsafe, and state 827# is shared across all uses (both in the user's code and in the Python 828# libraries), but that's fine for most programs and is easier for the 829# casual user than making them instantiate their own Random() instance. 830_inst = Random() 831seed = _inst.seed 832random = _inst.random 833uniform = _inst.uniform 834randint = _inst.randint 835choice = _inst.choice 836randrange = _inst.randrange 837sample = _inst.sample 838shuffle = _inst.shuffle 839normalvariate = _inst.normalvariate 840lognormvariate = _inst.lognormvariate 841cunifvariate = _inst.cunifvariate 842expovariate = _inst.expovariate 843vonmisesvariate = _inst.vonmisesvariate 844gammavariate = _inst.gammavariate 845stdgamma = _inst.stdgamma 846gauss = _inst.gauss 847betavariate = _inst.betavariate 848paretovariate = _inst.paretovariate 849weibullvariate = _inst.weibullvariate 850getstate = _inst.getstate 851setstate = _inst.setstate 852jumpahead = _inst.jumpahead 853whseed = _inst.whseed 854 855if __name__ == '__main__': 856 _test() 857