1"""Random variable generators. 2 3 bytes 4 ----- 5 uniform bytes (values between 0 and 255) 6 7 integers 8 -------- 9 uniform within range 10 11 sequences 12 --------- 13 pick random element 14 pick random sample 15 pick weighted random sample 16 generate random permutation 17 18 distributions on the real line: 19 ------------------------------ 20 uniform 21 triangular 22 normal (Gaussian) 23 lognormal 24 negative exponential 25 gamma 26 beta 27 pareto 28 Weibull 29 30 distributions on the circle (angles 0 to 2pi) 31 --------------------------------------------- 32 circular uniform 33 von Mises 34 35 discrete distributions 36 ---------------------- 37 binomial 38 39 40General notes on the underlying Mersenne Twister core generator: 41 42* The period is 2**19937-1. 43* It is one of the most extensively tested generators in existence. 44* The random() method is implemented in C, executes in a single Python step, 45 and is, therefore, threadsafe. 46 47""" 48 49# Translated by Guido van Rossum from C source provided by 50# Adrian Baddeley. Adapted by Raymond Hettinger for use with 51# the Mersenne Twister and os.urandom() core generators. 52 53from math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil 54from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin 55from math import tau as TWOPI, floor as _floor, isfinite as _isfinite 56from math import lgamma as _lgamma, fabs as _fabs, log2 as _log2 57from os import urandom as _urandom 58from _collections_abc import Sequence as _Sequence 59from operator import index as _index 60from itertools import accumulate as _accumulate, repeat as _repeat 61from bisect import bisect as _bisect 62import os as _os 63import _random 64 65__all__ = [ 66 "Random", 67 "SystemRandom", 68 "betavariate", 69 "binomialvariate", 70 "choice", 71 "choices", 72 "expovariate", 73 "gammavariate", 74 "gauss", 75 "getrandbits", 76 "getstate", 77 "lognormvariate", 78 "normalvariate", 79 "paretovariate", 80 "randbytes", 81 "randint", 82 "random", 83 "randrange", 84 "sample", 85 "seed", 86 "setstate", 87 "shuffle", 88 "triangular", 89 "uniform", 90 "vonmisesvariate", 91 "weibullvariate", 92] 93 94NV_MAGICCONST = 4 * _exp(-0.5) / _sqrt(2.0) 95LOG4 = _log(4.0) 96SG_MAGICCONST = 1.0 + _log(4.5) 97BPF = 53 # Number of bits in a float 98RECIP_BPF = 2 ** -BPF 99_ONE = 1 100_sha512 = None 101 102 103class Random(_random.Random): 104 """Random number generator base class used by bound module functions. 105 106 Used to instantiate instances of Random to get generators that don't 107 share state. 108 109 Class Random can also be subclassed if you want to use a different basic 110 generator of your own devising: in that case, override the following 111 methods: random(), seed(), getstate(), and setstate(). 112 Optionally, implement a getrandbits() method so that randrange() 113 can cover arbitrarily large ranges. 114 115 """ 116 117 VERSION = 3 # used by getstate/setstate 118 119 def __init__(self, x=None): 120 """Initialize an instance. 121 122 Optional argument x controls seeding, as for Random.seed(). 123 """ 124 125 self.seed(x) 126 self.gauss_next = None 127 128 def seed(self, a=None, version=2): 129 """Initialize internal state from a seed. 130 131 The only supported seed types are None, int, float, 132 str, bytes, and bytearray. 133 134 None or no argument seeds from current time or from an operating 135 system specific randomness source if available. 136 137 If *a* is an int, all bits are used. 138 139 For version 2 (the default), all of the bits are used if *a* is a str, 140 bytes, or bytearray. For version 1 (provided for reproducing random 141 sequences from older versions of Python), the algorithm for str and 142 bytes generates a narrower range of seeds. 143 144 """ 145 146 if version == 1 and isinstance(a, (str, bytes)): 147 a = a.decode('latin-1') if isinstance(a, bytes) else a 148 x = ord(a[0]) << 7 if a else 0 149 for c in map(ord, a): 150 x = ((1000003 * x) ^ c) & 0xFFFFFFFFFFFFFFFF 151 x ^= len(a) 152 a = -2 if x == -1 else x 153 154 elif version == 2 and isinstance(a, (str, bytes, bytearray)): 155 global _sha512 156 if _sha512 is None: 157 try: 158 # hashlib is pretty heavy to load, try lean internal 159 # module first 160 from _sha2 import sha512 as _sha512 161 except ImportError: 162 # fallback to official implementation 163 from hashlib import sha512 as _sha512 164 165 if isinstance(a, str): 166 a = a.encode() 167 a = int.from_bytes(a + _sha512(a).digest()) 168 169 elif not isinstance(a, (type(None), int, float, str, bytes, bytearray)): 170 raise TypeError('The only supported seed types are:\n' 171 'None, int, float, str, bytes, and bytearray.') 172 173 super().seed(a) 174 self.gauss_next = None 175 176 def getstate(self): 177 """Return internal state; can be passed to setstate() later.""" 178 return self.VERSION, super().getstate(), self.gauss_next 179 180 def setstate(self, state): 181 """Restore internal state from object returned by getstate().""" 182 version = state[0] 183 if version == 3: 184 version, internalstate, self.gauss_next = state 185 super().setstate(internalstate) 186 elif version == 2: 187 version, internalstate, self.gauss_next = state 188 # In version 2, the state was saved as signed ints, which causes 189 # inconsistencies between 32/64-bit systems. The state is 190 # really unsigned 32-bit ints, so we convert negative ints from 191 # version 2 to positive longs for version 3. 192 try: 193 internalstate = tuple(x % (2 ** 32) for x in internalstate) 194 except ValueError as e: 195 raise TypeError from e 196 super().setstate(internalstate) 197 else: 198 raise ValueError("state with version %s passed to " 199 "Random.setstate() of version %s" % 200 (version, self.VERSION)) 201 202 203 ## ------------------------------------------------------- 204 ## ---- Methods below this point do not need to be overridden or extended 205 ## ---- when subclassing for the purpose of using a different core generator. 206 207 208 ## -------------------- pickle support ------------------- 209 210 # Issue 17489: Since __reduce__ was defined to fix #759889 this is no 211 # longer called; we leave it here because it has been here since random was 212 # rewritten back in 2001 and why risk breaking something. 213 def __getstate__(self): # for pickle 214 return self.getstate() 215 216 def __setstate__(self, state): # for pickle 217 self.setstate(state) 218 219 def __reduce__(self): 220 return self.__class__, (), self.getstate() 221 222 223 ## ---- internal support method for evenly distributed integers ---- 224 225 def __init_subclass__(cls, /, **kwargs): 226 """Control how subclasses generate random integers. 227 228 The algorithm a subclass can use depends on the random() and/or 229 getrandbits() implementation available to it and determines 230 whether it can generate random integers from arbitrarily large 231 ranges. 232 """ 233 234 for c in cls.__mro__: 235 if '_randbelow' in c.__dict__: 236 # just inherit it 237 break 238 if 'getrandbits' in c.__dict__: 239 cls._randbelow = cls._randbelow_with_getrandbits 240 break 241 if 'random' in c.__dict__: 242 cls._randbelow = cls._randbelow_without_getrandbits 243 break 244 245 def _randbelow_with_getrandbits(self, n): 246 "Return a random int in the range [0,n). Defined for n > 0." 247 248 getrandbits = self.getrandbits 249 k = n.bit_length() 250 r = getrandbits(k) # 0 <= r < 2**k 251 while r >= n: 252 r = getrandbits(k) 253 return r 254 255 def _randbelow_without_getrandbits(self, n, maxsize=1<<BPF): 256 """Return a random int in the range [0,n). Defined for n > 0. 257 258 The implementation does not use getrandbits, but only random. 259 """ 260 261 random = self.random 262 if n >= maxsize: 263 from warnings import warn 264 warn("Underlying random() generator does not supply \n" 265 "enough bits to choose from a population range this large.\n" 266 "To remove the range limitation, add a getrandbits() method.") 267 return _floor(random() * n) 268 rem = maxsize % n 269 limit = (maxsize - rem) / maxsize # int(limit * maxsize) % n == 0 270 r = random() 271 while r >= limit: 272 r = random() 273 return _floor(r * maxsize) % n 274 275 _randbelow = _randbelow_with_getrandbits 276 277 278 ## -------------------------------------------------------- 279 ## ---- Methods below this point generate custom distributions 280 ## ---- based on the methods defined above. They do not 281 ## ---- directly touch the underlying generator and only 282 ## ---- access randomness through the methods: random(), 283 ## ---- getrandbits(), or _randbelow(). 284 285 286 ## -------------------- bytes methods --------------------- 287 288 def randbytes(self, n): 289 """Generate n random bytes.""" 290 return self.getrandbits(n * 8).to_bytes(n, 'little') 291 292 293 ## -------------------- integer methods ------------------- 294 295 def randrange(self, start, stop=None, step=_ONE): 296 """Choose a random item from range(stop) or range(start, stop[, step]). 297 298 Roughly equivalent to ``choice(range(start, stop, step))`` but 299 supports arbitrarily large ranges and is optimized for common cases. 300 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 = _index(start) 306 if stop is None: 307 # We don't check for "step != 1" because it hasn't been 308 # type checked and converted to an integer yet. 309 if step is not _ONE: 310 raise TypeError("Missing a non-None stop argument") 311 if istart > 0: 312 return self._randbelow(istart) 313 raise ValueError("empty range for randrange()") 314 315 # Stop argument supplied. 316 istop = _index(stop) 317 width = istop - istart 318 istep = _index(step) 319 # Fast path. 320 if istep == 1: 321 if width > 0: 322 return istart + self._randbelow(width) 323 raise ValueError(f"empty range in randrange({start}, {stop})") 324 325 # Non-unit step argument supplied. 326 if istep > 0: 327 n = (width + istep - 1) // istep 328 elif istep < 0: 329 n = (width + istep + 1) // istep 330 else: 331 raise ValueError("zero step for randrange()") 332 if n <= 0: 333 raise ValueError(f"empty range in randrange({start}, {stop}, {step})") 334 return istart + istep * self._randbelow(n) 335 336 def randint(self, a, b): 337 """Return random integer in range [a, b], including both end points. 338 """ 339 340 return self.randrange(a, b+1) 341 342 343 ## -------------------- sequence methods ------------------- 344 345 def choice(self, seq): 346 """Choose a random element from a non-empty sequence.""" 347 348 # As an accommodation for NumPy, we don't use "if not seq" 349 # because bool(numpy.array()) raises a ValueError. 350 if not len(seq): 351 raise IndexError('Cannot choose from an empty sequence') 352 return seq[self._randbelow(len(seq))] 353 354 def shuffle(self, x): 355 """Shuffle list x in place, and return None.""" 356 357 randbelow = self._randbelow 358 for i in reversed(range(1, len(x))): 359 # pick an element in x[:i+1] with which to exchange x[i] 360 j = randbelow(i + 1) 361 x[i], x[j] = x[j], x[i] 362 363 def sample(self, population, k, *, counts=None): 364 """Chooses k unique random elements from a population sequence. 365 366 Returns a new list containing elements from the population while 367 leaving the original population unchanged. The resulting list is 368 in selection order so that all sub-slices will also be valid random 369 samples. This allows raffle winners (the sample) to be partitioned 370 into grand prize and second place winners (the subslices). 371 372 Members of the population need not be hashable or unique. If the 373 population contains repeats, then each occurrence is a possible 374 selection in the sample. 375 376 Repeated elements can be specified one at a time or with the optional 377 counts parameter. For example: 378 379 sample(['red', 'blue'], counts=[4, 2], k=5) 380 381 is equivalent to: 382 383 sample(['red', 'red', 'red', 'red', 'blue', 'blue'], k=5) 384 385 To choose a sample from a range of integers, use range() for the 386 population argument. This is especially fast and space efficient 387 for sampling from a large population: 388 389 sample(range(10000000), 60) 390 391 """ 392 393 # Sampling without replacement entails tracking either potential 394 # selections (the pool) in a list or previous selections in a set. 395 396 # When the number of selections is small compared to the 397 # population, then tracking selections is efficient, requiring 398 # only a small set and an occasional reselection. For 399 # a larger number of selections, the pool tracking method is 400 # preferred since the list takes less space than the 401 # set and it doesn't suffer from frequent reselections. 402 403 # The number of calls to _randbelow() is kept at or near k, the 404 # theoretical minimum. This is important because running time 405 # is dominated by _randbelow() and because it extracts the 406 # least entropy from the underlying random number generators. 407 408 # Memory requirements are kept to the smaller of a k-length 409 # set or an n-length list. 410 411 # There are other sampling algorithms that do not require 412 # auxiliary memory, but they were rejected because they made 413 # too many calls to _randbelow(), making them slower and 414 # causing them to eat more entropy than necessary. 415 416 if not isinstance(population, _Sequence): 417 raise TypeError("Population must be a sequence. " 418 "For dicts or sets, use sorted(d).") 419 n = len(population) 420 if counts is not None: 421 cum_counts = list(_accumulate(counts)) 422 if len(cum_counts) != n: 423 raise ValueError('The number of counts does not match the population') 424 total = cum_counts.pop() 425 if not isinstance(total, int): 426 raise TypeError('Counts must be integers') 427 if total <= 0: 428 raise ValueError('Total of counts must be greater than zero') 429 selections = self.sample(range(total), k=k) 430 bisect = _bisect 431 return [population[bisect(cum_counts, s)] for s in selections] 432 randbelow = self._randbelow 433 if not 0 <= k <= n: 434 raise ValueError("Sample larger than population or is negative") 435 result = [None] * k 436 setsize = 21 # size of a small set minus size of an empty list 437 if k > 5: 438 setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets 439 if n <= setsize: 440 # An n-length list is smaller than a k-length set. 441 # Invariant: non-selected at pool[0 : n-i] 442 pool = list(population) 443 for i in range(k): 444 j = randbelow(n - i) 445 result[i] = pool[j] 446 pool[j] = pool[n - i - 1] # move non-selected item into vacancy 447 else: 448 selected = set() 449 selected_add = selected.add 450 for i in range(k): 451 j = randbelow(n) 452 while j in selected: 453 j = randbelow(n) 454 selected_add(j) 455 result[i] = population[j] 456 return result 457 458 def choices(self, population, weights=None, *, cum_weights=None, k=1): 459 """Return a k sized list of population elements chosen with replacement. 460 461 If the relative weights or cumulative weights are not specified, 462 the selections are made with equal probability. 463 464 """ 465 random = self.random 466 n = len(population) 467 if cum_weights is None: 468 if weights is None: 469 floor = _floor 470 n += 0.0 # convert to float for a small speed improvement 471 return [population[floor(random() * n)] for i in _repeat(None, k)] 472 try: 473 cum_weights = list(_accumulate(weights)) 474 except TypeError: 475 if not isinstance(weights, int): 476 raise 477 k = weights 478 raise TypeError( 479 f'The number of choices must be a keyword argument: {k=}' 480 ) from None 481 elif weights is not None: 482 raise TypeError('Cannot specify both weights and cumulative weights') 483 if len(cum_weights) != n: 484 raise ValueError('The number of weights does not match the population') 485 total = cum_weights[-1] + 0.0 # convert to float 486 if total <= 0.0: 487 raise ValueError('Total of weights must be greater than zero') 488 if not _isfinite(total): 489 raise ValueError('Total of weights must be finite') 490 bisect = _bisect 491 hi = n - 1 492 return [population[bisect(cum_weights, random() * total, 0, hi)] 493 for i in _repeat(None, k)] 494 495 496 ## -------------------- real-valued distributions ------------------- 497 498 def uniform(self, a, b): 499 """Get a random number in the range [a, b) or [a, b] depending on rounding. 500 501 The mean (expected value) and variance of the random variable are: 502 503 E[X] = (a + b) / 2 504 Var[X] = (b - a) ** 2 / 12 505 506 """ 507 return a + (b - a) * self.random() 508 509 def triangular(self, low=0.0, high=1.0, mode=None): 510 """Triangular distribution. 511 512 Continuous distribution bounded by given lower and upper limits, 513 and having a given mode value in-between. 514 515 http://en.wikipedia.org/wiki/Triangular_distribution 516 517 The mean (expected value) and variance of the random variable are: 518 519 E[X] = (low + high + mode) / 3 520 Var[X] = (low**2 + high**2 + mode**2 - low*high - low*mode - high*mode) / 18 521 522 """ 523 u = self.random() 524 try: 525 c = 0.5 if mode is None else (mode - low) / (high - low) 526 except ZeroDivisionError: 527 return low 528 if u > c: 529 u = 1.0 - u 530 c = 1.0 - c 531 low, high = high, low 532 return low + (high - low) * _sqrt(u * c) 533 534 def normalvariate(self, mu=0.0, sigma=1.0): 535 """Normal distribution. 536 537 mu is the mean, and sigma is the standard deviation. 538 539 """ 540 # Uses Kinderman and Monahan method. Reference: Kinderman, 541 # A.J. and Monahan, J.F., "Computer generation of random 542 # variables using the ratio of uniform deviates", ACM Trans 543 # Math Software, 3, (1977), pp257-260. 544 545 random = self.random 546 while True: 547 u1 = random() 548 u2 = 1.0 - random() 549 z = NV_MAGICCONST * (u1 - 0.5) / u2 550 zz = z * z / 4.0 551 if zz <= -_log(u2): 552 break 553 return mu + z * sigma 554 555 def gauss(self, mu=0.0, sigma=1.0): 556 """Gaussian distribution. 557 558 mu is the mean, and sigma is the standard deviation. This is 559 slightly faster than the normalvariate() function. 560 561 Not thread-safe without a lock around calls. 562 563 """ 564 # When x and y are two variables from [0, 1), uniformly 565 # distributed, then 566 # 567 # cos(2*pi*x)*sqrt(-2*log(1-y)) 568 # sin(2*pi*x)*sqrt(-2*log(1-y)) 569 # 570 # are two *independent* variables with normal distribution 571 # (mu = 0, sigma = 1). 572 # (Lambert Meertens) 573 # (corrected version; bug discovered by Mike Miller, fixed by LM) 574 575 # Multithreading note: When two threads call this function 576 # simultaneously, it is possible that they will receive the 577 # same return value. The window is very small though. To 578 # avoid this, you have to use a lock around all calls. (I 579 # didn't want to slow this down in the serial case by using a 580 # lock here.) 581 582 random = self.random 583 z = self.gauss_next 584 self.gauss_next = None 585 if z is None: 586 x2pi = random() * TWOPI 587 g2rad = _sqrt(-2.0 * _log(1.0 - random())) 588 z = _cos(x2pi) * g2rad 589 self.gauss_next = _sin(x2pi) * g2rad 590 591 return mu + z * sigma 592 593 def lognormvariate(self, mu, sigma): 594 """Log normal distribution. 595 596 If you take the natural logarithm of this distribution, you'll get a 597 normal distribution with mean mu and standard deviation sigma. 598 mu can have any value, and sigma must be greater than zero. 599 600 """ 601 return _exp(self.normalvariate(mu, sigma)) 602 603 def expovariate(self, lambd=1.0): 604 """Exponential distribution. 605 606 lambd is 1.0 divided by the desired mean. It should be 607 nonzero. (The parameter would be called "lambda", but that is 608 a reserved word in Python.) Returned values range from 0 to 609 positive infinity if lambd is positive, and from negative 610 infinity to 0 if lambd is negative. 611 612 The mean (expected value) and variance of the random variable are: 613 614 E[X] = 1 / lambd 615 Var[X] = 1 / lambd ** 2 616 617 """ 618 # we use 1-random() instead of random() to preclude the 619 # possibility of taking the log of zero. 620 621 return -_log(1.0 - self.random()) / lambd 622 623 def vonmisesvariate(self, mu, kappa): 624 """Circular data distribution. 625 626 mu is the mean angle, expressed in radians between 0 and 2*pi, and 627 kappa is the concentration parameter, which must be greater than or 628 equal to zero. If kappa is equal to zero, this distribution reduces 629 to a uniform random angle over the range 0 to 2*pi. 630 631 """ 632 # Based upon an algorithm published in: Fisher, N.I., 633 # "Statistical Analysis of Circular Data", Cambridge 634 # University Press, 1993. 635 636 # Thanks to Magnus Kessler for a correction to the 637 # implementation of step 4. 638 639 random = self.random 640 if kappa <= 1e-6: 641 return TWOPI * random() 642 643 s = 0.5 / kappa 644 r = s + _sqrt(1.0 + s * s) 645 646 while True: 647 u1 = random() 648 z = _cos(_pi * u1) 649 650 d = z / (r + z) 651 u2 = random() 652 if u2 < 1.0 - d * d or u2 <= (1.0 - d) * _exp(d): 653 break 654 655 q = 1.0 / r 656 f = (q + z) / (1.0 + q * z) 657 u3 = random() 658 if u3 > 0.5: 659 theta = (mu + _acos(f)) % TWOPI 660 else: 661 theta = (mu - _acos(f)) % TWOPI 662 663 return theta 664 665 def gammavariate(self, alpha, beta): 666 """Gamma distribution. Not the gamma function! 667 668 Conditions on the parameters are alpha > 0 and beta > 0. 669 670 The probability distribution function is: 671 672 x ** (alpha - 1) * math.exp(-x / beta) 673 pdf(x) = -------------------------------------- 674 math.gamma(alpha) * beta ** alpha 675 676 The mean (expected value) and variance of the random variable are: 677 678 E[X] = alpha * beta 679 Var[X] = alpha * beta ** 2 680 681 """ 682 683 # Warning: a few older sources define the gamma distribution in terms 684 # of alpha > -1.0 685 if alpha <= 0.0 or beta <= 0.0: 686 raise ValueError('gammavariate: alpha and beta must be > 0.0') 687 688 random = self.random 689 if alpha > 1.0: 690 691 # Uses R.C.H. Cheng, "The generation of Gamma 692 # variables with non-integral shape parameters", 693 # Applied Statistics, (1977), 26, No. 1, p71-74 694 695 ainv = _sqrt(2.0 * alpha - 1.0) 696 bbb = alpha - LOG4 697 ccc = alpha + ainv 698 699 while True: 700 u1 = random() 701 if not 1e-7 < u1 < 0.9999999: 702 continue 703 u2 = 1.0 - random() 704 v = _log(u1 / (1.0 - u1)) / ainv 705 x = alpha * _exp(v) 706 z = u1 * u1 * u2 707 r = bbb + ccc * v - x 708 if r + SG_MAGICCONST - 4.5 * z >= 0.0 or r >= _log(z): 709 return x * beta 710 711 elif alpha == 1.0: 712 # expovariate(1/beta) 713 return -_log(1.0 - random()) * beta 714 715 else: 716 # alpha is between 0 and 1 (exclusive) 717 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 718 while True: 719 u = random() 720 b = (_e + alpha) / _e 721 p = b * u 722 if p <= 1.0: 723 x = p ** (1.0 / alpha) 724 else: 725 x = -_log((b - p) / alpha) 726 u1 = random() 727 if p > 1.0: 728 if u1 <= x ** (alpha - 1.0): 729 break 730 elif u1 <= _exp(-x): 731 break 732 return x * beta 733 734 def betavariate(self, alpha, beta): 735 """Beta distribution. 736 737 Conditions on the parameters are alpha > 0 and beta > 0. 738 Returned values range between 0 and 1. 739 740 The mean (expected value) and variance of the random variable are: 741 742 E[X] = alpha / (alpha + beta) 743 Var[X] = alpha * beta / ((alpha + beta)**2 * (alpha + beta + 1)) 744 745 """ 746 ## See 747 ## http://mail.python.org/pipermail/python-bugs-list/2001-January/003752.html 748 ## for Ivan Frohne's insightful analysis of why the original implementation: 749 ## 750 ## def betavariate(self, alpha, beta): 751 ## # Discrete Event Simulation in C, pp 87-88. 752 ## 753 ## y = self.expovariate(alpha) 754 ## z = self.expovariate(1.0/beta) 755 ## return z/(y+z) 756 ## 757 ## was dead wrong, and how it probably got that way. 758 759 # This version due to Janne Sinkkonen, and matches all the std 760 # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 761 y = self.gammavariate(alpha, 1.0) 762 if y: 763 return y / (y + self.gammavariate(beta, 1.0)) 764 return 0.0 765 766 def paretovariate(self, alpha): 767 """Pareto distribution. alpha is the shape parameter.""" 768 # Jain, pg. 495 769 770 u = 1.0 - self.random() 771 return u ** (-1.0 / alpha) 772 773 def weibullvariate(self, alpha, beta): 774 """Weibull distribution. 775 776 alpha is the scale parameter and beta is the shape parameter. 777 778 """ 779 # Jain, pg. 499; bug fix courtesy Bill Arms 780 781 u = 1.0 - self.random() 782 return alpha * (-_log(u)) ** (1.0 / beta) 783 784 785 ## -------------------- discrete distributions --------------------- 786 787 def binomialvariate(self, n=1, p=0.5): 788 """Binomial random variable. 789 790 Gives the number of successes for *n* independent trials 791 with the probability of success in each trial being *p*: 792 793 sum(random() < p for i in range(n)) 794 795 Returns an integer in the range: 0 <= X <= n 796 797 The mean (expected value) and variance of the random variable are: 798 799 E[X] = n * p 800 Var[x] = n * p * (1 - p) 801 802 """ 803 # Error check inputs and handle edge cases 804 if n < 0: 805 raise ValueError("n must be non-negative") 806 if p <= 0.0 or p >= 1.0: 807 if p == 0.0: 808 return 0 809 if p == 1.0: 810 return n 811 raise ValueError("p must be in the range 0.0 <= p <= 1.0") 812 813 random = self.random 814 815 # Fast path for a common case 816 if n == 1: 817 return _index(random() < p) 818 819 # Exploit symmetry to establish: p <= 0.5 820 if p > 0.5: 821 return n - self.binomialvariate(n, 1.0 - p) 822 823 if n * p < 10.0: 824 # BG: Geometric method by Devroye with running time of O(np). 825 # https://dl.acm.org/doi/pdf/10.1145/42372.42381 826 x = y = 0 827 c = _log2(1.0 - p) 828 if not c: 829 return x 830 while True: 831 y += _floor(_log2(random()) / c) + 1 832 if y > n: 833 return x 834 x += 1 835 836 # BTRS: Transformed rejection with squeeze method by Wolfgang Hörmann 837 # https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.47.8407&rep=rep1&type=pdf 838 assert n*p >= 10.0 and p <= 0.5 839 setup_complete = False 840 841 spq = _sqrt(n * p * (1.0 - p)) # Standard deviation of the distribution 842 b = 1.15 + 2.53 * spq 843 a = -0.0873 + 0.0248 * b + 0.01 * p 844 c = n * p + 0.5 845 vr = 0.92 - 4.2 / b 846 847 while True: 848 849 u = random() 850 u -= 0.5 851 us = 0.5 - _fabs(u) 852 k = _floor((2.0 * a / us + b) * u + c) 853 if k < 0 or k > n: 854 continue 855 856 # The early-out "squeeze" test substantially reduces 857 # the number of acceptance condition evaluations. 858 v = random() 859 if us >= 0.07 and v <= vr: 860 return k 861 862 # Acceptance-rejection test. 863 # Note, the original paper erroneously omits the call to log(v) 864 # when comparing to the log of the rescaled binomial distribution. 865 if not setup_complete: 866 alpha = (2.83 + 5.1 / b) * spq 867 lpq = _log(p / (1.0 - p)) 868 m = _floor((n + 1) * p) # Mode of the distribution 869 h = _lgamma(m + 1) + _lgamma(n - m + 1) 870 setup_complete = True # Only needs to be done once 871 v *= alpha / (a / (us * us) + b) 872 if _log(v) <= h - _lgamma(k + 1) - _lgamma(n - k + 1) + (k - m) * lpq: 873 return k 874 875 876## ------------------------------------------------------------------ 877## --------------- Operating System Random Source ------------------ 878 879 880class SystemRandom(Random): 881 """Alternate random number generator using sources provided 882 by the operating system (such as /dev/urandom on Unix or 883 CryptGenRandom on Windows). 884 885 Not available on all systems (see os.urandom() for details). 886 887 """ 888 889 def random(self): 890 """Get the next random number in the range 0.0 <= X < 1.0.""" 891 return (int.from_bytes(_urandom(7)) >> 3) * RECIP_BPF 892 893 def getrandbits(self, k): 894 """getrandbits(k) -> x. Generates an int with k random bits.""" 895 if k < 0: 896 raise ValueError('number of bits must be non-negative') 897 numbytes = (k + 7) // 8 # bits / 8 and rounded up 898 x = int.from_bytes(_urandom(numbytes)) 899 return x >> (numbytes * 8 - k) # trim excess bits 900 901 def randbytes(self, n): 902 """Generate n random bytes.""" 903 # os.urandom(n) fails with ValueError for n < 0 904 # and returns an empty bytes string for n == 0. 905 return _urandom(n) 906 907 def seed(self, *args, **kwds): 908 "Stub method. Not used for a system random number generator." 909 return None 910 911 def _notimplemented(self, *args, **kwds): 912 "Method should not be called for a system random number generator." 913 raise NotImplementedError('System entropy source does not have state.') 914 getstate = setstate = _notimplemented 915 916 917# ---------------------------------------------------------------------- 918# Create one instance, seeded from current time, and export its methods 919# as module-level functions. The functions share state across all uses 920# (both in the user's code and in the Python libraries), but that's fine 921# for most programs and is easier for the casual user than making them 922# instantiate their own Random() instance. 923 924_inst = Random() 925seed = _inst.seed 926random = _inst.random 927uniform = _inst.uniform 928triangular = _inst.triangular 929randint = _inst.randint 930choice = _inst.choice 931randrange = _inst.randrange 932sample = _inst.sample 933shuffle = _inst.shuffle 934choices = _inst.choices 935normalvariate = _inst.normalvariate 936lognormvariate = _inst.lognormvariate 937expovariate = _inst.expovariate 938vonmisesvariate = _inst.vonmisesvariate 939gammavariate = _inst.gammavariate 940gauss = _inst.gauss 941betavariate = _inst.betavariate 942binomialvariate = _inst.binomialvariate 943paretovariate = _inst.paretovariate 944weibullvariate = _inst.weibullvariate 945getstate = _inst.getstate 946setstate = _inst.setstate 947getrandbits = _inst.getrandbits 948randbytes = _inst.randbytes 949 950 951## ------------------------------------------------------ 952## ----------------- test program ----------------------- 953 954def _test_generator(n, func, args): 955 from statistics import stdev, fmean as mean 956 from time import perf_counter 957 958 t0 = perf_counter() 959 data = [func(*args) for i in _repeat(None, n)] 960 t1 = perf_counter() 961 962 xbar = mean(data) 963 sigma = stdev(data, xbar) 964 low = min(data) 965 high = max(data) 966 967 print(f'{t1 - t0:.3f} sec, {n} times {func.__name__}{args!r}') 968 print('avg %g, stddev %g, min %g, max %g\n' % (xbar, sigma, low, high)) 969 970 971def _test(N=10_000): 972 _test_generator(N, random, ()) 973 _test_generator(N, normalvariate, (0.0, 1.0)) 974 _test_generator(N, lognormvariate, (0.0, 1.0)) 975 _test_generator(N, vonmisesvariate, (0.0, 1.0)) 976 _test_generator(N, binomialvariate, (15, 0.60)) 977 _test_generator(N, binomialvariate, (100, 0.75)) 978 _test_generator(N, gammavariate, (0.01, 1.0)) 979 _test_generator(N, gammavariate, (0.1, 1.0)) 980 _test_generator(N, gammavariate, (0.1, 2.0)) 981 _test_generator(N, gammavariate, (0.5, 1.0)) 982 _test_generator(N, gammavariate, (0.9, 1.0)) 983 _test_generator(N, gammavariate, (1.0, 1.0)) 984 _test_generator(N, gammavariate, (2.0, 1.0)) 985 _test_generator(N, gammavariate, (20.0, 1.0)) 986 _test_generator(N, gammavariate, (200.0, 1.0)) 987 _test_generator(N, gauss, (0.0, 1.0)) 988 _test_generator(N, betavariate, (3.0, 3.0)) 989 _test_generator(N, triangular, (0.0, 1.0, 1.0 / 3.0)) 990 991 992## ------------------------------------------------------ 993## ------------------ fork support --------------------- 994 995if hasattr(_os, "fork"): 996 _os.register_at_fork(after_in_child=_inst.seed) 997 998 999# ------------------------------------------------------ 1000# -------------- command-line interface ---------------- 1001 1002 1003def _parse_args(arg_list: list[str] | None): 1004 import argparse 1005 parser = argparse.ArgumentParser( 1006 formatter_class=argparse.RawTextHelpFormatter) 1007 group = parser.add_mutually_exclusive_group() 1008 group.add_argument( 1009 "-c", "--choice", nargs="+", 1010 help="print a random choice") 1011 group.add_argument( 1012 "-i", "--integer", type=int, metavar="N", 1013 help="print a random integer between 1 and N inclusive") 1014 group.add_argument( 1015 "-f", "--float", type=float, metavar="N", 1016 help="print a random floating-point number between 0 and N inclusive") 1017 group.add_argument( 1018 "--test", type=int, const=10_000, nargs="?", 1019 help=argparse.SUPPRESS) 1020 parser.add_argument("input", nargs="*", 1021 help="""\ 1022if no options given, output depends on the input 1023 string or multiple: same as --choice 1024 integer: same as --integer 1025 float: same as --float""") 1026 args = parser.parse_args(arg_list) 1027 return args, parser.format_help() 1028 1029 1030def main(arg_list: list[str] | None = None) -> int | str: 1031 args, help_text = _parse_args(arg_list) 1032 1033 # Explicit arguments 1034 if args.choice: 1035 return choice(args.choice) 1036 1037 if args.integer is not None: 1038 return randint(1, args.integer) 1039 1040 if args.float is not None: 1041 return uniform(0, args.float) 1042 1043 if args.test: 1044 _test(args.test) 1045 return "" 1046 1047 # No explicit argument, select based on input 1048 if len(args.input) == 1: 1049 val = args.input[0] 1050 try: 1051 # Is it an integer? 1052 val = int(val) 1053 return randint(1, val) 1054 except ValueError: 1055 try: 1056 # Is it a float? 1057 val = float(val) 1058 return uniform(0, val) 1059 except ValueError: 1060 # Split in case of space-separated string: "a b c" 1061 return choice(val.split()) 1062 1063 if len(args.input) >= 2: 1064 return choice(args.input) 1065 1066 return help_text 1067 1068 1069if __name__ == '__main__': 1070 print(main()) 1071