blob: 66433baa8cc18f2781c082c15ac97184dd14693a [file] [log] [blame]
Guido van Rossume7b146f2000-02-04 15:28:42 +00001"""Random variable generators.
Guido van Rossumff03b1a1994-03-09 12:55:02 +00002
Raymond Hettingeref19bad2020-06-25 17:03:50 -07003 bytes
4 -----
5 uniform bytes (values between 0 and 255)
6
Tim Petersd7b5e882001-01-25 03:36:26 +00007 integers
8 --------
9 uniform within range
10
11 sequences
12 ---------
13 pick random element
Raymond Hettingerf24eb352002-11-12 17:41:57 +000014 pick random sample
Raymond Hettingere8f1e002016-09-06 17:15:29 -070015 pick weighted random sample
Tim Petersd7b5e882001-01-25 03:36:26 +000016 generate random permutation
17
Guido van Rossume7b146f2000-02-04 15:28:42 +000018 distributions on the real line:
19 ------------------------------
Tim Petersd7b5e882001-01-25 03:36:26 +000020 uniform
Christian Heimesfe337bf2008-03-23 21:54:12 +000021 triangular
Guido van Rossume7b146f2000-02-04 15:28:42 +000022 normal (Gaussian)
23 lognormal
24 negative exponential
25 gamma
26 beta
Raymond Hettinger40f62172002-12-29 23:03:38 +000027 pareto
28 Weibull
Guido van Rossumff03b1a1994-03-09 12:55:02 +000029
Guido van Rossume7b146f2000-02-04 15:28:42 +000030 distributions on the circle (angles 0 to 2pi)
31 ---------------------------------------------
32 circular uniform
33 von Mises
34
Raymond Hettinger40f62172002-12-29 23:03:38 +000035General notes on the underlying Mersenne Twister core generator:
Guido van Rossume7b146f2000-02-04 15:28:42 +000036
Raymond Hettinger40f62172002-12-29 23:03:38 +000037* The period is 2**19937-1.
Thomas Wouters0e3f5912006-08-11 14:57:12 +000038* It is one of the most extensively tested generators in existence.
Thomas Wouters0e3f5912006-08-11 14:57:12 +000039* The random() method is implemented in C, executes in a single Python step,
40 and is, therefore, threadsafe.
Tim Peterse360d952001-01-26 10:00:39 +000041
Guido van Rossume7b146f2000-02-04 15:28:42 +000042"""
Guido van Rossumd03e1191998-05-29 17:51:31 +000043
Raymond Hettingeref19bad2020-06-25 17:03:50 -070044# Translated by Guido van Rossum from C source provided by
45# Adrian Baddeley. Adapted by Raymond Hettinger for use with
46# the Mersenne Twister and os.urandom() core generators.
47
Raymond Hettinger2f726e92003-10-05 09:09:15 +000048from warnings import warn as _warn
Raymond Hettinger91e27c22005-08-19 01:36:35 +000049from math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil
Raymond Hettinger26a1ad12020-06-22 19:38:59 -070050from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
Ram Rachumb0dfc752020-09-29 04:32:10 +030051from math import tau as TWOPI, floor as _floor, isfinite as _isfinite
Raymond Hettingerc1c43ca2004-09-05 00:00:42 +000052from os import urandom as _urandom
Christian Heimesf1dc3ee2013-10-13 02:04:20 +020053from _collections_abc import Set as _Set, Sequence as _Sequence
Raymond Hettinger53822032019-02-16 13:30:51 -080054from itertools import accumulate as _accumulate, repeat as _repeat
Raymond Hettingercfd31f02019-02-13 02:04:17 -080055from bisect import bisect as _bisect
Antoine Pitrou346cbd32017-05-27 17:50:54 +020056import os as _os
Raymond Hettingeref19bad2020-06-25 17:03:50 -070057import _random
Guido van Rossumff03b1a1994-03-09 12:55:02 +000058
Christian Heimesd9145962019-04-10 22:18:02 +020059try:
60 # hashlib is pretty heavy to load, try lean internal module first
61 from _sha512 import sha512 as _sha512
62except ImportError:
63 # fallback to official implementation
64 from hashlib import sha512 as _sha512
65
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -070066__all__ = [
67 "Random",
68 "SystemRandom",
69 "betavariate",
70 "choice",
71 "choices",
72 "expovariate",
73 "gammavariate",
74 "gauss",
75 "getrandbits",
76 "getstate",
77 "lognormvariate",
78 "normalvariate",
79 "paretovariate",
80 "randint",
81 "random",
82 "randrange",
83 "sample",
84 "seed",
85 "setstate",
86 "shuffle",
87 "triangular",
88 "uniform",
89 "vonmisesvariate",
90 "weibullvariate",
91]
Tim Petersd7b5e882001-01-25 03:36:26 +000092
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -070093NV_MAGICCONST = 4 * _exp(-0.5) / _sqrt(2.0)
Tim Petersd7b5e882001-01-25 03:36:26 +000094LOG4 = _log(4.0)
Tim Petersd7b5e882001-01-25 03:36:26 +000095SG_MAGICCONST = 1.0 + _log(4.5)
Raymond Hettinger2f726e92003-10-05 09:09:15 +000096BPF = 53 # Number of bits in a float
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -070097RECIP_BPF = 2 ** -BPF
Tim Petersd7b5e882001-01-25 03:36:26 +000098
Raymond Hettinger356a4592004-08-30 06:14:31 +000099
Raymond Hettinger145a4a02003-01-07 10:25:55 +0000100class Random(_random.Random):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000101 """Random number generator base class used by bound module functions.
102
103 Used to instantiate instances of Random to get generators that don't
Raymond Hettinger28de64f2008-01-13 23:40:30 +0000104 share state.
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000105
106 Class Random can also be subclassed if you want to use a different basic
107 generator of your own devising: in that case, override the following
Raymond Hettinger28de64f2008-01-13 23:40:30 +0000108 methods: random(), seed(), getstate(), and setstate().
Benjamin Petersond18de0e2008-07-31 20:21:46 +0000109 Optionally, implement a getrandbits() method so that randrange()
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000110 can cover arbitrarily large ranges.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000111
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000112 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000113
Christian Heimescbf3b5c2007-12-03 21:02:03 +0000114 VERSION = 3 # used by getstate/setstate
Tim Petersd7b5e882001-01-25 03:36:26 +0000115
116 def __init__(self, x=None):
117 """Initialize an instance.
118
119 Optional argument x controls seeding, as for Random.seed().
120 """
121
122 self.seed(x)
Raymond Hettinger40f62172002-12-29 23:03:38 +0000123 self.gauss_next = None
Tim Petersd7b5e882001-01-25 03:36:26 +0000124
Raymond Hettingerf763a722010-09-07 00:38:15 +0000125 def seed(self, a=None, version=2):
Raymond Hettingerd0cdeaa2019-08-22 09:19:36 -0700126 """Initialize internal state from a seed.
127
128 The only supported seed types are None, int, float,
129 str, bytes, and bytearray.
Tim Petersd7b5e882001-01-25 03:36:26 +0000130
Raymond Hettinger23f12412004-09-13 22:23:21 +0000131 None or no argument seeds from current time or from an operating
132 system specific randomness source if available.
Tim Peters0de88fc2001-02-01 04:59:18 +0000133
Raymond Hettinger183cd1f2010-09-08 18:48:21 +0000134 If *a* is an int, all bits are used.
Raymond Hettingerf763a722010-09-07 00:38:15 +0000135
Raymond Hettinger16eb8272016-09-04 11:17:28 -0700136 For version 2 (the default), all of the bits are used if *a* is a str,
137 bytes, or bytearray. For version 1 (provided for reproducing random
138 sequences from older versions of Python), the algorithm for str and
139 bytes generates a narrower range of seeds.
140
Tim Petersd7b5e882001-01-25 03:36:26 +0000141 """
142
Raymond Hettingerc7bab7c2016-08-31 15:01:08 -0700143 if version == 1 and isinstance(a, (str, bytes)):
Raymond Hettinger132a7d72017-09-17 09:04:30 -0700144 a = a.decode('latin-1') if isinstance(a, bytes) else a
Raymond Hettingerc7bab7c2016-08-31 15:01:08 -0700145 x = ord(a[0]) << 7 if a else 0
Raymond Hettinger132a7d72017-09-17 09:04:30 -0700146 for c in map(ord, a):
147 x = ((1000003 * x) ^ c) & 0xFFFFFFFFFFFFFFFF
Raymond Hettingerc7bab7c2016-08-31 15:01:08 -0700148 x ^= len(a)
149 a = -2 if x == -1 else x
150
Raymond Hettingerd0cdeaa2019-08-22 09:19:36 -0700151 elif version == 2 and isinstance(a, (str, bytes, bytearray)):
Raymond Hettinger2f9cc7a2016-08-31 23:00:32 -0700152 if isinstance(a, str):
153 a = a.encode()
154 a += _sha512(a).digest()
155 a = int.from_bytes(a, 'big')
Raymond Hettingerf763a722010-09-07 00:38:15 +0000156
Raymond Hettingerd0cdeaa2019-08-22 09:19:36 -0700157 elif not isinstance(a, (type(None), int, float, str, bytes, bytearray)):
158 _warn('Seeding based on hashing is deprecated\n'
159 'since Python 3.9 and will be removed in a subsequent '
160 'version. The only \n'
161 'supported seed types are: None, '
162 'int, float, str, bytes, and bytearray.',
163 DeprecationWarning, 2)
164
Guido van Rossumcd16bf62007-06-13 18:07:49 +0000165 super().seed(a)
Tim Peters46c04e12002-05-05 20:40:00 +0000166 self.gauss_next = None
167
Tim Peterscd804102001-01-25 20:25:57 +0000168 def getstate(self):
169 """Return internal state; can be passed to setstate() later."""
Guido van Rossumcd16bf62007-06-13 18:07:49 +0000170 return self.VERSION, super().getstate(), self.gauss_next
Tim Peterscd804102001-01-25 20:25:57 +0000171
172 def setstate(self, state):
173 """Restore internal state from object returned by getstate()."""
174 version = state[0]
Christian Heimescbf3b5c2007-12-03 21:02:03 +0000175 if version == 3:
Raymond Hettinger40f62172002-12-29 23:03:38 +0000176 version, internalstate, self.gauss_next = state
Guido van Rossumcd16bf62007-06-13 18:07:49 +0000177 super().setstate(internalstate)
Christian Heimescbf3b5c2007-12-03 21:02:03 +0000178 elif version == 2:
179 version, internalstate, self.gauss_next = state
180 # In version 2, the state was saved as signed ints, which causes
181 # inconsistencies between 32/64-bit systems. The state is
182 # really unsigned 32-bit ints, so we convert negative ints from
183 # version 2 to positive longs for version 3.
184 try:
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700185 internalstate = tuple(x % (2 ** 32) for x in internalstate)
Christian Heimescbf3b5c2007-12-03 21:02:03 +0000186 except ValueError as e:
187 raise TypeError from e
Raymond Hettinger183cd1f2010-09-08 18:48:21 +0000188 super().setstate(internalstate)
Tim Peterscd804102001-01-25 20:25:57 +0000189 else:
190 raise ValueError("state with version %s passed to "
191 "Random.setstate() of version %s" %
192 (version, self.VERSION))
193
Tim Peterscd804102001-01-25 20:25:57 +0000194
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700195 ## -------------------------------------------------------
196 ## ---- Methods below this point do not need to be overridden or extended
197 ## ---- when subclassing for the purpose of using a different core generator.
Victor Stinner2d875772020-04-29 18:49:00 +0200198
Victor Stinner2d875772020-04-29 18:49:00 +0200199
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700200 ## -------------------- pickle support -------------------
Tim Peterscd804102001-01-25 20:25:57 +0000201
R David Murrayd9ebf4d2013-04-02 13:10:52 -0400202 # Issue 17489: Since __reduce__ was defined to fix #759889 this is no
203 # longer called; we leave it here because it has been here since random was
204 # rewritten back in 2001 and why risk breaking something.
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700205 def __getstate__(self): # for pickle
Tim Peterscd804102001-01-25 20:25:57 +0000206 return self.getstate()
207
208 def __setstate__(self, state): # for pickle
209 self.setstate(state)
210
Raymond Hettinger5f078ff2003-06-24 20:29:04 +0000211 def __reduce__(self):
212 return self.__class__, (), self.getstate()
213
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700214
215 ## ---- internal support method for evenly distributed integers ----
216
217 def __init_subclass__(cls, /, **kwargs):
218 """Control how subclasses generate random integers.
219
220 The algorithm a subclass can use depends on the random() and/or
221 getrandbits() implementation available to it and determines
222 whether it can generate random integers from arbitrarily large
223 ranges.
224 """
225
226 for c in cls.__mro__:
227 if '_randbelow' in c.__dict__:
228 # just inherit it
229 break
230 if 'getrandbits' in c.__dict__:
231 cls._randbelow = cls._randbelow_with_getrandbits
232 break
233 if 'random' in c.__dict__:
234 cls._randbelow = cls._randbelow_without_getrandbits
235 break
236
237 def _randbelow_with_getrandbits(self, n):
238 "Return a random int in the range [0,n). Returns 0 if n==0."
239
240 if not n:
241 return 0
242 getrandbits = self.getrandbits
243 k = n.bit_length() # don't use (n-1) here because n can be 1
244 r = getrandbits(k) # 0 <= r < 2**k
245 while r >= n:
246 r = getrandbits(k)
247 return r
248
249 def _randbelow_without_getrandbits(self, n, maxsize=1<<BPF):
250 """Return a random int in the range [0,n). Returns 0 if n==0.
251
252 The implementation does not use getrandbits, but only random.
253 """
254
255 random = self.random
256 if n >= maxsize:
257 _warn("Underlying random() generator does not supply \n"
258 "enough bits to choose from a population range this large.\n"
259 "To remove the range limitation, add a getrandbits() method.")
260 return _floor(random() * n)
261 if n == 0:
262 return 0
263 rem = maxsize % n
264 limit = (maxsize - rem) / maxsize # int(limit * maxsize) % n == 0
265 r = random()
266 while r >= limit:
267 r = random()
268 return _floor(r * maxsize) % n
269
270 _randbelow = _randbelow_with_getrandbits
271
272
273 ## --------------------------------------------------------
274 ## ---- Methods below this point generate custom distributions
275 ## ---- based on the methods defined above. They do not
276 ## ---- directly touch the underlying generator and only
277 ## ---- access randomness through the methods: random(),
278 ## ---- getrandbits(), or _randbelow().
279
280
281 ## -------------------- bytes methods ---------------------
282
283 def randbytes(self, n):
284 """Generate n random bytes."""
285 return self.getrandbits(n * 8).to_bytes(n, 'little')
286
287
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700288 ## -------------------- integer methods -------------------
Tim Peterscd804102001-01-25 20:25:57 +0000289
Raymond Hettinger26a1ad12020-06-22 19:38:59 -0700290 def randrange(self, start, stop=None, step=1):
Tim Petersd7b5e882001-01-25 03:36:26 +0000291 """Choose a random item from range(start, stop[, step]).
292
293 This fixes the problem with randint() which includes the
294 endpoint; in Python this is usually not what you want.
Raymond Hettinger3051cc32010-09-07 00:48:40 +0000295
Tim Petersd7b5e882001-01-25 03:36:26 +0000296 """
297
298 # This code is a bit messy to make it fast for the
Tim Peters9146f272002-08-16 03:41:39 +0000299 # common case while still doing adequate error checking.
Raymond Hettinger26a1ad12020-06-22 19:38:59 -0700300 istart = int(start)
Tim Petersd7b5e882001-01-25 03:36:26 +0000301 if istart != start:
Collin Winterce36ad82007-08-30 01:19:48 +0000302 raise ValueError("non-integer arg 1 for randrange()")
Raymond Hettinger3051cc32010-09-07 00:48:40 +0000303 if stop is None:
Tim Petersd7b5e882001-01-25 03:36:26 +0000304 if istart > 0:
Raymond Hettinger05156612010-09-07 04:44:52 +0000305 return self._randbelow(istart)
Collin Winterce36ad82007-08-30 01:19:48 +0000306 raise ValueError("empty range for randrange()")
Tim Peters9146f272002-08-16 03:41:39 +0000307
308 # stop argument supplied.
Raymond Hettinger26a1ad12020-06-22 19:38:59 -0700309 istop = int(stop)
Tim Petersd7b5e882001-01-25 03:36:26 +0000310 if istop != stop:
Collin Winterce36ad82007-08-30 01:19:48 +0000311 raise ValueError("non-integer stop for randrange()")
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000312 width = istop - istart
313 if step == 1 and width > 0:
Raymond Hettingerc3246972010-09-07 09:32:57 +0000314 return istart + self._randbelow(width)
Tim Petersd7b5e882001-01-25 03:36:26 +0000315 if step == 1:
Kumar Akshay2433a2a2019-01-22 00:49:59 +0530316 raise ValueError("empty range for randrange() (%d, %d, %d)" % (istart, istop, width))
Tim Peters9146f272002-08-16 03:41:39 +0000317
318 # Non-unit step argument supplied.
Raymond Hettinger26a1ad12020-06-22 19:38:59 -0700319 istep = int(step)
Tim Petersd7b5e882001-01-25 03:36:26 +0000320 if istep != step:
Collin Winterce36ad82007-08-30 01:19:48 +0000321 raise ValueError("non-integer step for randrange()")
Tim Petersd7b5e882001-01-25 03:36:26 +0000322 if istep > 0:
Raymond Hettingerffdb8bb2004-09-27 15:29:05 +0000323 n = (width + istep - 1) // istep
Tim Petersd7b5e882001-01-25 03:36:26 +0000324 elif istep < 0:
Raymond Hettingerffdb8bb2004-09-27 15:29:05 +0000325 n = (width + istep + 1) // istep
Tim Petersd7b5e882001-01-25 03:36:26 +0000326 else:
Collin Winterce36ad82007-08-30 01:19:48 +0000327 raise ValueError("zero step for randrange()")
Tim Petersd7b5e882001-01-25 03:36:26 +0000328
329 if n <= 0:
Collin Winterce36ad82007-08-30 01:19:48 +0000330 raise ValueError("empty range for randrange()")
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000331
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700332 return istart + istep * self._randbelow(n)
Tim Petersd7b5e882001-01-25 03:36:26 +0000333
334 def randint(self, a, b):
Tim Peterscd804102001-01-25 20:25:57 +0000335 """Return random integer in range [a, b], including both end points.
Tim Petersd7b5e882001-01-25 03:36:26 +0000336 """
337
338 return self.randrange(a, b+1)
339
Wolfgang Maierba3a87a2018-04-17 17:16:17 +0200340
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700341 ## -------------------- sequence methods -------------------
Tim Peterscd804102001-01-25 20:25:57 +0000342
Tim Petersd7b5e882001-01-25 03:36:26 +0000343 def choice(self, seq):
344 """Choose a random element from a non-empty sequence."""
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700345 # raises IndexError if seq is empty
346 return seq[self._randbelow(len(seq))]
Tim Petersd7b5e882001-01-25 03:36:26 +0000347
Raymond Hettinger8fe47c32013-10-05 21:48:21 -0700348 def shuffle(self, x, random=None):
Antoine Pitrou5e394332012-11-04 02:10:33 +0100349 """Shuffle list x in place, and return None.
Tim Petersd7b5e882001-01-25 03:36:26 +0000350
Antoine Pitrou5e394332012-11-04 02:10:33 +0100351 Optional argument random is a 0-argument function returning a
352 random float in [0.0, 1.0); if it is the default None, the
353 standard random.random will be used.
Senthil Kumaranf8ce51a2013-09-11 22:54:31 -0700354
Tim Petersd7b5e882001-01-25 03:36:26 +0000355 """
356
Raymond Hettinger8fe47c32013-10-05 21:48:21 -0700357 if random is None:
358 randbelow = self._randbelow
359 for i in reversed(range(1, len(x))):
360 # pick an element in x[:i+1] with which to exchange x[i]
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700361 j = randbelow(i + 1)
Raymond Hettinger8fe47c32013-10-05 21:48:21 -0700362 x[i], x[j] = x[j], x[i]
363 else:
Raymond Hettinger190fac92020-05-02 16:45:32 -0700364 _warn('The *random* parameter to shuffle() has been deprecated\n'
365 'since Python 3.9 and will be removed in a subsequent '
366 'version.',
367 DeprecationWarning, 2)
Raymond Hettinger26a1ad12020-06-22 19:38:59 -0700368 floor = _floor
Raymond Hettinger8fe47c32013-10-05 21:48:21 -0700369 for i in reversed(range(1, len(x))):
370 # pick an element in x[:i+1] with which to exchange x[i]
Raymond Hettinger26a1ad12020-06-22 19:38:59 -0700371 j = floor(random() * (i + 1))
Raymond Hettinger8fe47c32013-10-05 21:48:21 -0700372 x[i], x[j] = x[j], x[i]
Tim Petersd7b5e882001-01-25 03:36:26 +0000373
Raymond Hettinger81a5fc32020-05-08 07:53:15 -0700374 def sample(self, population, k, *, counts=None):
Raymond Hettinger1acde192008-01-14 01:00:53 +0000375 """Chooses k unique random elements from a population sequence or set.
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000376
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000377 Returns a new list containing elements from the population while
378 leaving the original population unchanged. The resulting list is
379 in selection order so that all sub-slices will also be valid random
380 samples. This allows raffle winners (the sample) to be partitioned
381 into grand prize and second place winners (the subslices).
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000382
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000383 Members of the population need not be hashable or unique. If the
384 population contains repeats, then each occurrence is a possible
385 selection in the sample.
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000386
Raymond Hettinger81a5fc32020-05-08 07:53:15 -0700387 Repeated elements can be specified one at a time or with the optional
388 counts parameter. For example:
389
390 sample(['red', 'blue'], counts=[4, 2], k=5)
391
392 is equivalent to:
393
394 sample(['red', 'red', 'red', 'red', 'blue', 'blue'], k=5)
395
396 To choose a sample from a range of integers, use range() for the
397 population argument. This is especially fast and space efficient
398 for sampling from a large population:
399
400 sample(range(10000000), 60)
401
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000402 """
403
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000404 # Sampling without replacement entails tracking either potential
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000405 # selections (the pool) in a list or previous selections in a set.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000406
Jeremy Hylton2b55d352004-02-23 17:27:57 +0000407 # When the number of selections is small compared to the
408 # population, then tracking selections is efficient, requiring
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000409 # only a small set and an occasional reselection. For
Jeremy Hylton2b55d352004-02-23 17:27:57 +0000410 # a larger number of selections, the pool tracking method is
411 # preferred since the list takes less space than the
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000412 # set and it doesn't suffer from frequent reselections.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000413
Raymond Hettinger7fc633f2018-12-04 00:13:38 -0800414 # The number of calls to _randbelow() is kept at or near k, the
415 # theoretical minimum. This is important because running time
416 # is dominated by _randbelow() and because it extracts the
417 # least entropy from the underlying random number generators.
418
419 # Memory requirements are kept to the smaller of a k-length
420 # set or an n-length list.
421
422 # There are other sampling algorithms that do not require
423 # auxiliary memory, but they were rejected because they made
424 # too many calls to _randbelow(), making them slower and
425 # causing them to eat more entropy than necessary.
426
Raymond Hettinger57d1a882011-02-23 00:46:28 +0000427 if not isinstance(population, _Sequence):
masklinn1e27b572020-12-19 05:33:36 +0100428 if isinstance(population, _Set):
429 _warn('Sampling from a set deprecated\n'
430 'since Python 3.9 and will be removed in a subsequent version.',
431 DeprecationWarning, 2)
432 population = tuple(population)
433 else:
434 raise TypeError("Population must be a sequence. For dicts or sets, use sorted(d).")
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000435 n = len(population)
Raymond Hettinger81a5fc32020-05-08 07:53:15 -0700436 if counts is not None:
437 cum_counts = list(_accumulate(counts))
438 if len(cum_counts) != n:
439 raise ValueError('The number of counts does not match the population')
440 total = cum_counts.pop()
441 if not isinstance(total, int):
442 raise TypeError('Counts must be integers')
443 if total <= 0:
444 raise ValueError('Total of counts must be greater than zero')
445 selections = sample(range(total), k=k)
446 bisect = _bisect
447 return [population[bisect(cum_counts, s)] for s in selections]
448 randbelow = self._randbelow
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000449 if not 0 <= k <= n:
Raymond Hettingerbf871262016-11-21 14:34:33 -0800450 raise ValueError("Sample larger than population or is negative")
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000451 result = [None] * k
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000452 setsize = 21 # size of a small set minus size of an empty list
453 if k > 5:
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700454 setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets
Raymond Hettinger1acde192008-01-14 01:00:53 +0000455 if n <= setsize:
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700456 # An n-length list is smaller than a k-length set.
457 # Invariant: non-selected at pool[0 : n-i]
Raymond Hettinger311f4192002-11-18 09:01:24 +0000458 pool = list(population)
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700459 for i in range(k):
460 j = randbelow(n - i)
Raymond Hettinger311f4192002-11-18 09:01:24 +0000461 result[i] = pool[j]
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700462 pool[j] = pool[n - i - 1] # move non-selected item into vacancy
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000463 else:
Raymond Hettinger1acde192008-01-14 01:00:53 +0000464 selected = set()
465 selected_add = selected.add
466 for i in range(k):
Raymond Hettinger05a505f2010-09-07 19:19:33 +0000467 j = randbelow(n)
Raymond Hettinger1acde192008-01-14 01:00:53 +0000468 while j in selected:
Raymond Hettinger05a505f2010-09-07 19:19:33 +0000469 j = randbelow(n)
Raymond Hettinger1acde192008-01-14 01:00:53 +0000470 selected_add(j)
471 result[i] = population[j]
Raymond Hettinger311f4192002-11-18 09:01:24 +0000472 return result
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000473
Raymond Hettinger9016f282016-09-26 21:45:57 -0700474 def choices(self, population, weights=None, *, cum_weights=None, k=1):
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700475 """Return a k sized list of population elements chosen with replacement.
476
477 If the relative weights or cumulative weights are not specified,
478 the selections are made with equal probability.
479
480 """
Raymond Hettinger30d00e52016-10-29 16:55:36 -0700481 random = self.random
Raymond Hettingere69cd162018-07-04 15:28:20 -0700482 n = len(population)
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700483 if cum_weights is None:
484 if weights is None:
Raymond Hettinger26a1ad12020-06-22 19:38:59 -0700485 floor = _floor
Raymond Hettinger0a18e052018-11-09 02:39:50 -0800486 n += 0.0 # convert to float for a small speed improvement
Raymond Hettinger26a1ad12020-06-22 19:38:59 -0700487 return [population[floor(random() * n)] for i in _repeat(None, k)]
Raymond Hettingercfd31f02019-02-13 02:04:17 -0800488 cum_weights = list(_accumulate(weights))
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700489 elif weights is not None:
Raymond Hettinger24e42392016-11-13 00:42:56 -0500490 raise TypeError('Cannot specify both weights and cumulative weights')
Raymond Hettingere69cd162018-07-04 15:28:20 -0700491 if len(cum_weights) != n:
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700492 raise ValueError('The number of weights does not match the population')
Raymond Hettinger0a18e052018-11-09 02:39:50 -0800493 total = cum_weights[-1] + 0.0 # convert to float
Raymond Hettinger041d8b42019-11-23 02:22:13 -0800494 if total <= 0.0:
495 raise ValueError('Total of weights must be greater than zero')
Ram Rachumb0dfc752020-09-29 04:32:10 +0300496 if not _isfinite(total):
497 raise ValueError('Total of weights must be finite')
Raymond Hettinger041d8b42019-11-23 02:22:13 -0800498 bisect = _bisect
Raymond Hettingere69cd162018-07-04 15:28:20 -0700499 hi = n - 1
Raymond Hettingerddf71712018-06-27 01:08:31 -0700500 return [population[bisect(cum_weights, random() * total, 0, hi)]
Raymond Hettinger53822032019-02-16 13:30:51 -0800501 for i in _repeat(None, k)]
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700502
Tim Peterscd804102001-01-25 20:25:57 +0000503
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700504 ## -------------------- real-valued distributions -------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000505
506 def uniform(self, a, b):
Raymond Hettingerbe40db02009-06-11 23:12:14 +0000507 "Get a random number in the range [a, b) or [a, b] depending on rounding."
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700508 return a + (b - a) * self.random()
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000509
Christian Heimesfe337bf2008-03-23 21:54:12 +0000510 def triangular(self, low=0.0, high=1.0, mode=None):
511 """Triangular distribution.
512
513 Continuous distribution bounded by given lower and upper limits,
514 and having a given mode value in-between.
515
516 http://en.wikipedia.org/wiki/Triangular_distribution
517
518 """
519 u = self.random()
Raymond Hettinger978c6ab2014-05-25 17:25:27 -0700520 try:
521 c = 0.5 if mode is None else (mode - low) / (high - low)
522 except ZeroDivisionError:
523 return low
Christian Heimesfe337bf2008-03-23 21:54:12 +0000524 if u > c:
525 u = 1.0 - u
526 c = 1.0 - c
527 low, high = high, low
Raymond Hettingerf5ea83f2017-09-04 16:51:06 -0700528 return low + (high - low) * _sqrt(u * c)
Christian Heimesfe337bf2008-03-23 21:54:12 +0000529
Tim Petersd7b5e882001-01-25 03:36:26 +0000530 def normalvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000531 """Normal distribution.
532
533 mu is the mean, and sigma is the standard deviation.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000534
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000535 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000536 # Uses Kinderman and Monahan method. Reference: Kinderman,
537 # A.J. and Monahan, J.F., "Computer generation of random
538 # variables using the ratio of uniform deviates", ACM Trans
539 # Math Software, 3, (1977), pp257-260.
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000540
Tim Petersd7b5e882001-01-25 03:36:26 +0000541 random = self.random
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700542 while True:
Tim Peters0c9886d2001-01-15 01:18:21 +0000543 u1 = random()
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000544 u2 = 1.0 - random()
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700545 z = NV_MAGICCONST * (u1 - 0.5) / u2
546 zz = z * z / 4.0
Tim Petersd7b5e882001-01-25 03:36:26 +0000547 if zz <= -_log(u2):
548 break
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700549 return mu + z * sigma
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000550
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700551 def gauss(self, mu, sigma):
552 """Gaussian distribution.
553
554 mu is the mean, and sigma is the standard deviation. This is
555 slightly faster than the normalvariate() function.
556
557 Not thread-safe without a lock around calls.
558
559 """
560 # When x and y are two variables from [0, 1), uniformly
561 # distributed, then
562 #
563 # cos(2*pi*x)*sqrt(-2*log(1-y))
564 # sin(2*pi*x)*sqrt(-2*log(1-y))
565 #
566 # are two *independent* variables with normal distribution
567 # (mu = 0, sigma = 1).
568 # (Lambert Meertens)
569 # (corrected version; bug discovered by Mike Miller, fixed by LM)
570
571 # Multithreading note: When two threads call this function
572 # simultaneously, it is possible that they will receive the
573 # same return value. The window is very small though. To
574 # avoid this, you have to use a lock around all calls. (I
575 # didn't want to slow this down in the serial case by using a
576 # lock here.)
577
578 random = self.random
579 z = self.gauss_next
580 self.gauss_next = None
581 if z is None:
582 x2pi = random() * TWOPI
583 g2rad = _sqrt(-2.0 * _log(1.0 - random()))
584 z = _cos(x2pi) * g2rad
585 self.gauss_next = _sin(x2pi) * g2rad
586
587 return mu + z * sigma
Tim Petersd7b5e882001-01-25 03:36:26 +0000588
589 def lognormvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000590 """Log normal distribution.
591
592 If you take the natural logarithm of this distribution, you'll get a
593 normal distribution with mean mu and standard deviation sigma.
594 mu can have any value, and sigma must be greater than zero.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000595
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000596 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000597 return _exp(self.normalvariate(mu, sigma))
598
Tim Petersd7b5e882001-01-25 03:36:26 +0000599 def expovariate(self, lambd):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000600 """Exponential distribution.
601
Mark Dickinson2f947362009-01-07 17:54:07 +0000602 lambd is 1.0 divided by the desired mean. It should be
603 nonzero. (The parameter would be called "lambda", but that is
604 a reserved word in Python.) Returned values range from 0 to
605 positive infinity if lambd is positive, and from negative
606 infinity to 0 if lambd is negative.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000607
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000608 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000609 # lambd: rate lambd = 1/mean
610 # ('lambda' is a Python reserved word)
611
Raymond Hettinger5279fb92011-06-25 11:30:53 +0200612 # we use 1-random() instead of random() to preclude the
613 # possibility of taking the log of zero.
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700614 return -_log(1.0 - self.random()) / lambd
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000615
Tim Petersd7b5e882001-01-25 03:36:26 +0000616 def vonmisesvariate(self, mu, kappa):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000617 """Circular data distribution.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000618
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000619 mu is the mean angle, expressed in radians between 0 and 2*pi, and
620 kappa is the concentration parameter, which must be greater than or
621 equal to zero. If kappa is equal to zero, this distribution reduces
622 to a uniform random angle over the range 0 to 2*pi.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000623
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000624 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000625 # Based upon an algorithm published in: Fisher, N.I.,
626 # "Statistical Analysis of Circular Data", Cambridge
627 # University Press, 1993.
628
629 # Thanks to Magnus Kessler for a correction to the
630 # implementation of step 4.
631
632 random = self.random
633 if kappa <= 1e-6:
634 return TWOPI * random()
635
Serhiy Storchaka6c22b1d2013-02-10 19:28:56 +0200636 s = 0.5 / kappa
637 r = s + _sqrt(1.0 + s * s)
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000638
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700639 while True:
Tim Peters0c9886d2001-01-15 01:18:21 +0000640 u1 = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000641 z = _cos(_pi * u1)
Tim Petersd7b5e882001-01-25 03:36:26 +0000642
Serhiy Storchaka6c22b1d2013-02-10 19:28:56 +0200643 d = z / (r + z)
Tim Petersd7b5e882001-01-25 03:36:26 +0000644 u2 = random()
Serhiy Storchaka6c22b1d2013-02-10 19:28:56 +0200645 if u2 < 1.0 - d * d or u2 <= (1.0 - d) * _exp(d):
Tim Peters0c9886d2001-01-15 01:18:21 +0000646 break
Tim Petersd7b5e882001-01-25 03:36:26 +0000647
Serhiy Storchaka6c22b1d2013-02-10 19:28:56 +0200648 q = 1.0 / r
649 f = (q + z) / (1.0 + q * z)
Tim Petersd7b5e882001-01-25 03:36:26 +0000650 u3 = random()
651 if u3 > 0.5:
Mark Dickinsonbe5f9192013-02-10 14:16:10 +0000652 theta = (mu + _acos(f)) % TWOPI
Tim Petersd7b5e882001-01-25 03:36:26 +0000653 else:
Mark Dickinsonbe5f9192013-02-10 14:16:10 +0000654 theta = (mu - _acos(f)) % TWOPI
Tim Petersd7b5e882001-01-25 03:36:26 +0000655
656 return theta
657
Tim Petersd7b5e882001-01-25 03:36:26 +0000658 def gammavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000659 """Gamma distribution. Not the gamma function!
660
661 Conditions on the parameters are alpha > 0 and beta > 0.
662
Raymond Hettingera8e4d6e2011-03-22 15:55:51 -0700663 The probability distribution function is:
664
665 x ** (alpha - 1) * math.exp(-x / beta)
666 pdf(x) = --------------------------------------
667 math.gamma(alpha) * beta ** alpha
668
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000669 """
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000670 # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
Tim Peters8ac14952002-05-23 15:15:30 +0000671
Guido van Rossum570764d2002-05-14 14:08:12 +0000672 # Warning: a few older sources define the gamma distribution in terms
673 # of alpha > -1.0
674 if alpha <= 0.0 or beta <= 0.0:
Collin Winterce36ad82007-08-30 01:19:48 +0000675 raise ValueError('gammavariate: alpha and beta must be > 0.0')
Tim Peters8ac14952002-05-23 15:15:30 +0000676
Tim Petersd7b5e882001-01-25 03:36:26 +0000677 random = self.random
Tim Petersd7b5e882001-01-25 03:36:26 +0000678 if alpha > 1.0:
679
680 # Uses R.C.H. Cheng, "The generation of Gamma
681 # variables with non-integral shape parameters",
682 # Applied Statistics, (1977), 26, No. 1, p71-74
683
Raymond Hettingerca6cdc22002-05-13 23:40:14 +0000684 ainv = _sqrt(2.0 * alpha - 1.0)
685 bbb = alpha - LOG4
686 ccc = alpha + ainv
Tim Peters8ac14952002-05-23 15:15:30 +0000687
Raymond Hettinger6a613f92020-08-02 12:03:32 -0700688 while True:
Tim Petersd7b5e882001-01-25 03:36:26 +0000689 u1 = random()
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700690 if not 1e-7 < u1 < 0.9999999:
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000691 continue
692 u2 = 1.0 - random()
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700693 v = _log(u1 / (1.0 - u1)) / ainv
694 x = alpha * _exp(v)
695 z = u1 * u1 * u2
696 r = bbb + ccc * v - x
697 if r + SG_MAGICCONST - 4.5 * z >= 0.0 or r >= _log(z):
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000698 return x * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000699
700 elif alpha == 1.0:
leodema9f396b62017-06-04 07:41:41 +0100701 # expovariate(1/beta)
leodema63d15222018-12-24 07:54:25 +0100702 return -_log(1.0 - random()) * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000703
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700704 else:
705 # alpha is between 0 and 1 (exclusive)
Tim Petersd7b5e882001-01-25 03:36:26 +0000706 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700707 while True:
Tim Petersd7b5e882001-01-25 03:36:26 +0000708 u = random()
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700709 b = (_e + alpha) / _e
710 p = b * u
Tim Petersd7b5e882001-01-25 03:36:26 +0000711 if p <= 1.0:
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700712 x = p ** (1.0 / alpha)
Tim Petersd7b5e882001-01-25 03:36:26 +0000713 else:
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700714 x = -_log((b - p) / alpha)
Tim Petersd7b5e882001-01-25 03:36:26 +0000715 u1 = random()
Raymond Hettinger42406e62005-04-30 09:02:51 +0000716 if p > 1.0:
717 if u1 <= x ** (alpha - 1.0):
718 break
719 elif u1 <= _exp(-x):
Tim Petersd7b5e882001-01-25 03:36:26 +0000720 break
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000721 return x * beta
722
Tim Petersd7b5e882001-01-25 03:36:26 +0000723 def betavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000724 """Beta distribution.
725
Thomas Woutersb2137042007-02-01 18:02:27 +0000726 Conditions on the parameters are alpha > 0 and beta > 0.
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000727 Returned values range between 0 and 1.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000728
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000729 """
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700730 ## See
731 ## http://mail.python.org/pipermail/python-bugs-list/2001-January/003752.html
732 ## for Ivan Frohne's insightful analysis of why the original implementation:
733 ##
734 ## def betavariate(self, alpha, beta):
735 ## # Discrete Event Simulation in C, pp 87-88.
736 ##
737 ## y = self.expovariate(alpha)
738 ## z = self.expovariate(1.0/beta)
739 ## return z/(y+z)
740 ##
741 ## was dead wrong, and how it probably got that way.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000742
Tim Peters85e2e472001-01-26 06:49:56 +0000743 # This version due to Janne Sinkkonen, and matches all the std
744 # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
Raymond Hettinger650c1c92016-06-25 05:36:42 +0300745 y = self.gammavariate(alpha, 1.0)
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700746 if y:
Raymond Hettinger650c1c92016-06-25 05:36:42 +0300747 return y / (y + self.gammavariate(beta, 1.0))
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700748 return 0.0
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000749
Tim Petersd7b5e882001-01-25 03:36:26 +0000750 def paretovariate(self, alpha):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000751 """Pareto distribution. alpha is the shape parameter."""
Tim Petersd7b5e882001-01-25 03:36:26 +0000752 # Jain, pg. 495
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000753
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000754 u = 1.0 - self.random()
Raymond Hettinger5c327092020-08-01 01:18:26 -0700755 return u ** (-1.0 / alpha)
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000756
Tim Petersd7b5e882001-01-25 03:36:26 +0000757 def weibullvariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000758 """Weibull distribution.
759
760 alpha is the scale parameter and beta is the shape parameter.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000761
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000762 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000763 # Jain, pg. 499; bug fix courtesy Bill Arms
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000764
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000765 u = 1.0 - self.random()
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700766 return alpha * (-_log(u)) ** (1.0 / beta)
767
Guido van Rossum6c395ba1999-08-18 13:53:28 +0000768
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700769## ------------------------------------------------------------------
Raymond Hettinger23f12412004-09-13 22:23:21 +0000770## --------------- Operating System Random Source ------------------
Raymond Hettinger356a4592004-08-30 06:14:31 +0000771
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700772
Raymond Hettinger23f12412004-09-13 22:23:21 +0000773class SystemRandom(Random):
774 """Alternate random number generator using sources provided
775 by the operating system (such as /dev/urandom on Unix or
776 CryptGenRandom on Windows).
Raymond Hettinger356a4592004-08-30 06:14:31 +0000777
778 Not available on all systems (see os.urandom() for details).
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700779
Raymond Hettinger356a4592004-08-30 06:14:31 +0000780 """
781
782 def random(self):
783 """Get the next random number in the range [0.0, 1.0)."""
Raymond Hettinger183cd1f2010-09-08 18:48:21 +0000784 return (int.from_bytes(_urandom(7), 'big') >> 3) * RECIP_BPF
Raymond Hettinger356a4592004-08-30 06:14:31 +0000785
786 def getrandbits(self, k):
Serhiy Storchaka95949422013-08-27 19:40:23 +0300787 """getrandbits(k) -> x. Generates an int with k random bits."""
Antoine Pitrou75a33782020-04-17 19:32:14 +0200788 if k < 0:
789 raise ValueError('number of bits must be non-negative')
Raymond Hettinger63b17672010-09-08 19:27:59 +0000790 numbytes = (k + 7) // 8 # bits / 8 and rounded up
791 x = int.from_bytes(_urandom(numbytes), 'big')
792 return x >> (numbytes * 8 - k) # trim excess bits
Raymond Hettinger356a4592004-08-30 06:14:31 +0000793
Victor Stinner9f5fe792020-04-17 19:05:35 +0200794 def randbytes(self, n):
795 """Generate n random bytes."""
796 # os.urandom(n) fails with ValueError for n < 0
797 # and returns an empty bytes string for n == 0.
798 return _urandom(n)
799
Raymond Hettinger28de64f2008-01-13 23:40:30 +0000800 def seed(self, *args, **kwds):
Raymond Hettinger23f12412004-09-13 22:23:21 +0000801 "Stub method. Not used for a system random number generator."
Raymond Hettinger356a4592004-08-30 06:14:31 +0000802 return None
Raymond Hettinger356a4592004-08-30 06:14:31 +0000803
804 def _notimplemented(self, *args, **kwds):
Raymond Hettinger23f12412004-09-13 22:23:21 +0000805 "Method should not be called for a system random number generator."
806 raise NotImplementedError('System entropy source does not have state.')
Raymond Hettinger356a4592004-08-30 06:14:31 +0000807 getstate = setstate = _notimplemented
808
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700809
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700810# ----------------------------------------------------------------------
811# Create one instance, seeded from current time, and export its methods
812# as module-level functions. The functions share state across all uses
813# (both in the user's code and in the Python libraries), but that's fine
814# for most programs and is easier for the casual user than making them
815# instantiate their own Random() instance.
816
817_inst = Random()
818seed = _inst.seed
819random = _inst.random
820uniform = _inst.uniform
821triangular = _inst.triangular
822randint = _inst.randint
823choice = _inst.choice
824randrange = _inst.randrange
825sample = _inst.sample
826shuffle = _inst.shuffle
827choices = _inst.choices
828normalvariate = _inst.normalvariate
829lognormvariate = _inst.lognormvariate
830expovariate = _inst.expovariate
831vonmisesvariate = _inst.vonmisesvariate
832gammavariate = _inst.gammavariate
833gauss = _inst.gauss
834betavariate = _inst.betavariate
835paretovariate = _inst.paretovariate
836weibullvariate = _inst.weibullvariate
837getstate = _inst.getstate
838setstate = _inst.setstate
839getrandbits = _inst.getrandbits
840randbytes = _inst.randbytes
841
842
843## ------------------------------------------------------
844## ----------------- test program -----------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000845
Raymond Hettinger62297132003-08-30 01:24:19 +0000846def _test_generator(n, func, args):
Raymond Hettinger26a1ad12020-06-22 19:38:59 -0700847 from statistics import stdev, fmean as mean
848 from time import perf_counter
849
850 t0 = perf_counter()
851 data = [func(*args) for i in range(n)]
852 t1 = perf_counter()
853
854 xbar = mean(data)
855 sigma = stdev(data, xbar)
856 low = min(data)
857 high = max(data)
858
859 print(f'{t1 - t0:.3f} sec, {n} times {func.__name__}')
860 print('avg %g, stddev %g, min %g, max %g\n' % (xbar, sigma, low, high))
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000861
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000862
863def _test(N=2000):
Raymond Hettinger62297132003-08-30 01:24:19 +0000864 _test_generator(N, random, ())
865 _test_generator(N, normalvariate, (0.0, 1.0))
866 _test_generator(N, lognormvariate, (0.0, 1.0))
867 _test_generator(N, vonmisesvariate, (0.0, 1.0))
868 _test_generator(N, gammavariate, (0.01, 1.0))
869 _test_generator(N, gammavariate, (0.1, 1.0))
870 _test_generator(N, gammavariate, (0.1, 2.0))
871 _test_generator(N, gammavariate, (0.5, 1.0))
872 _test_generator(N, gammavariate, (0.9, 1.0))
873 _test_generator(N, gammavariate, (1.0, 1.0))
874 _test_generator(N, gammavariate, (2.0, 1.0))
875 _test_generator(N, gammavariate, (20.0, 1.0))
876 _test_generator(N, gammavariate, (200.0, 1.0))
877 _test_generator(N, gauss, (0.0, 1.0))
878 _test_generator(N, betavariate, (3.0, 3.0))
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700879 _test_generator(N, triangular, (0.0, 1.0, 1.0 / 3.0))
Tim Peterscd804102001-01-25 20:25:57 +0000880
Raymond Hettinger40f62172002-12-29 23:03:38 +0000881
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700882## ------------------------------------------------------
883## ------------------ fork support ---------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000884
Antoine Pitrou346cbd32017-05-27 17:50:54 +0200885if hasattr(_os, "fork"):
Gregory P. Smith163468a2017-05-29 10:03:41 -0700886 _os.register_at_fork(after_in_child=_inst.seed)
Antoine Pitrou346cbd32017-05-27 17:50:54 +0200887
888
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000889if __name__ == '__main__':
Tim Petersd7b5e882001-01-25 03:36:26 +0000890 _test()