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