blob: 139e8a40bb2724b4eea784dcac30310403138a94 [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 isinstance(population, _Set):
Raymond Hettinger4fe00202020-04-19 00:36:42 -0700428 _warn('Sampling from a set deprecated\n'
429 'since Python 3.9 and will be removed in a subsequent version.',
430 DeprecationWarning, 2)
Raymond Hettinger1acde192008-01-14 01:00:53 +0000431 population = tuple(population)
Raymond Hettinger57d1a882011-02-23 00:46:28 +0000432 if not isinstance(population, _Sequence):
Raymond Hettinger4fe00202020-04-19 00:36:42 -0700433 raise TypeError("Population must be a sequence. For dicts or sets, use sorted(d).")
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000434 n = len(population)
Raymond Hettinger81a5fc32020-05-08 07:53:15 -0700435 if counts is not None:
436 cum_counts = list(_accumulate(counts))
437 if len(cum_counts) != n:
438 raise ValueError('The number of counts does not match the population')
439 total = cum_counts.pop()
440 if not isinstance(total, int):
441 raise TypeError('Counts must be integers')
442 if total <= 0:
443 raise ValueError('Total of counts must be greater than zero')
444 selections = sample(range(total), k=k)
445 bisect = _bisect
446 return [population[bisect(cum_counts, s)] for s in selections]
447 randbelow = self._randbelow
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000448 if not 0 <= k <= n:
Raymond Hettingerbf871262016-11-21 14:34:33 -0800449 raise ValueError("Sample larger than population or is negative")
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000450 result = [None] * k
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000451 setsize = 21 # size of a small set minus size of an empty list
452 if k > 5:
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700453 setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets
Raymond Hettinger1acde192008-01-14 01:00:53 +0000454 if n <= setsize:
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700455 # An n-length list is smaller than a k-length set.
456 # Invariant: non-selected at pool[0 : n-i]
Raymond Hettinger311f4192002-11-18 09:01:24 +0000457 pool = list(population)
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700458 for i in range(k):
459 j = randbelow(n - i)
Raymond Hettinger311f4192002-11-18 09:01:24 +0000460 result[i] = pool[j]
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700461 pool[j] = pool[n - i - 1] # move non-selected item into vacancy
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000462 else:
Raymond Hettinger1acde192008-01-14 01:00:53 +0000463 selected = set()
464 selected_add = selected.add
465 for i in range(k):
Raymond Hettinger05a505f2010-09-07 19:19:33 +0000466 j = randbelow(n)
Raymond Hettinger1acde192008-01-14 01:00:53 +0000467 while j in selected:
Raymond Hettinger05a505f2010-09-07 19:19:33 +0000468 j = randbelow(n)
Raymond Hettinger1acde192008-01-14 01:00:53 +0000469 selected_add(j)
470 result[i] = population[j]
Raymond Hettinger311f4192002-11-18 09:01:24 +0000471 return result
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000472
Raymond Hettinger9016f282016-09-26 21:45:57 -0700473 def choices(self, population, weights=None, *, cum_weights=None, k=1):
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700474 """Return a k sized list of population elements chosen with replacement.
475
476 If the relative weights or cumulative weights are not specified,
477 the selections are made with equal probability.
478
479 """
Raymond Hettinger30d00e52016-10-29 16:55:36 -0700480 random = self.random
Raymond Hettingere69cd162018-07-04 15:28:20 -0700481 n = len(population)
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700482 if cum_weights is None:
483 if weights is None:
Raymond Hettinger26a1ad12020-06-22 19:38:59 -0700484 floor = _floor
Raymond Hettinger0a18e052018-11-09 02:39:50 -0800485 n += 0.0 # convert to float for a small speed improvement
Raymond Hettinger26a1ad12020-06-22 19:38:59 -0700486 return [population[floor(random() * n)] for i in _repeat(None, k)]
Raymond Hettingercfd31f02019-02-13 02:04:17 -0800487 cum_weights = list(_accumulate(weights))
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700488 elif weights is not None:
Raymond Hettinger24e42392016-11-13 00:42:56 -0500489 raise TypeError('Cannot specify both weights and cumulative weights')
Raymond Hettingere69cd162018-07-04 15:28:20 -0700490 if len(cum_weights) != n:
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700491 raise ValueError('The number of weights does not match the population')
Raymond Hettinger0a18e052018-11-09 02:39:50 -0800492 total = cum_weights[-1] + 0.0 # convert to float
Raymond Hettinger041d8b42019-11-23 02:22:13 -0800493 if total <= 0.0:
494 raise ValueError('Total of weights must be greater than zero')
Ram Rachumb0dfc752020-09-29 04:32:10 +0300495 if not _isfinite(total):
496 raise ValueError('Total of weights must be finite')
Raymond Hettinger041d8b42019-11-23 02:22:13 -0800497 bisect = _bisect
Raymond Hettingere69cd162018-07-04 15:28:20 -0700498 hi = n - 1
Raymond Hettingerddf71712018-06-27 01:08:31 -0700499 return [population[bisect(cum_weights, random() * total, 0, hi)]
Raymond Hettinger53822032019-02-16 13:30:51 -0800500 for i in _repeat(None, k)]
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700501
Tim Peterscd804102001-01-25 20:25:57 +0000502
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700503 ## -------------------- real-valued distributions -------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000504
505 def uniform(self, a, b):
Raymond Hettingerbe40db02009-06-11 23:12:14 +0000506 "Get a random number in the range [a, b) or [a, b] depending on rounding."
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700507 return a + (b - a) * self.random()
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000508
Christian Heimesfe337bf2008-03-23 21:54:12 +0000509 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 """
518 u = self.random()
Raymond Hettinger978c6ab2014-05-25 17:25:27 -0700519 try:
520 c = 0.5 if mode is None else (mode - low) / (high - low)
521 except ZeroDivisionError:
522 return low
Christian Heimesfe337bf2008-03-23 21:54:12 +0000523 if u > c:
524 u = 1.0 - u
525 c = 1.0 - c
526 low, high = high, low
Raymond Hettingerf5ea83f2017-09-04 16:51:06 -0700527 return low + (high - low) * _sqrt(u * c)
Christian Heimesfe337bf2008-03-23 21:54:12 +0000528
Tim Petersd7b5e882001-01-25 03:36:26 +0000529 def normalvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000530 """Normal distribution.
531
532 mu is the mean, and sigma is the standard deviation.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000533
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000534 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000535 # Uses Kinderman and Monahan method. Reference: Kinderman,
536 # A.J. and Monahan, J.F., "Computer generation of random
537 # variables using the ratio of uniform deviates", ACM Trans
538 # Math Software, 3, (1977), pp257-260.
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000539
Tim Petersd7b5e882001-01-25 03:36:26 +0000540 random = self.random
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700541 while True:
Tim Peters0c9886d2001-01-15 01:18:21 +0000542 u1 = random()
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000543 u2 = 1.0 - random()
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700544 z = NV_MAGICCONST * (u1 - 0.5) / u2
545 zz = z * z / 4.0
Tim Petersd7b5e882001-01-25 03:36:26 +0000546 if zz <= -_log(u2):
547 break
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700548 return mu + z * sigma
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000549
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700550 def gauss(self, mu, sigma):
551 """Gaussian distribution.
552
553 mu is the mean, and sigma is the standard deviation. This is
554 slightly faster than the normalvariate() function.
555
556 Not thread-safe without a lock around calls.
557
558 """
559 # When x and y are two variables from [0, 1), uniformly
560 # distributed, then
561 #
562 # cos(2*pi*x)*sqrt(-2*log(1-y))
563 # sin(2*pi*x)*sqrt(-2*log(1-y))
564 #
565 # are two *independent* variables with normal distribution
566 # (mu = 0, sigma = 1).
567 # (Lambert Meertens)
568 # (corrected version; bug discovered by Mike Miller, fixed by LM)
569
570 # Multithreading note: When two threads call this function
571 # simultaneously, it is possible that they will receive the
572 # same return value. The window is very small though. To
573 # avoid this, you have to use a lock around all calls. (I
574 # didn't want to slow this down in the serial case by using a
575 # lock here.)
576
577 random = self.random
578 z = self.gauss_next
579 self.gauss_next = None
580 if z is None:
581 x2pi = random() * TWOPI
582 g2rad = _sqrt(-2.0 * _log(1.0 - random()))
583 z = _cos(x2pi) * g2rad
584 self.gauss_next = _sin(x2pi) * g2rad
585
586 return mu + z * sigma
Tim Petersd7b5e882001-01-25 03:36:26 +0000587
588 def lognormvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000589 """Log normal distribution.
590
591 If you take the natural logarithm of this distribution, you'll get a
592 normal distribution with mean mu and standard deviation sigma.
593 mu can have any value, and sigma must be greater than zero.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000594
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000595 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000596 return _exp(self.normalvariate(mu, sigma))
597
Tim Petersd7b5e882001-01-25 03:36:26 +0000598 def expovariate(self, lambd):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000599 """Exponential distribution.
600
Mark Dickinson2f947362009-01-07 17:54:07 +0000601 lambd is 1.0 divided by the desired mean. It should be
602 nonzero. (The parameter would be called "lambda", but that is
603 a reserved word in Python.) Returned values range from 0 to
604 positive infinity if lambd is positive, and from negative
605 infinity to 0 if lambd is negative.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000606
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000607 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000608 # lambd: rate lambd = 1/mean
609 # ('lambda' is a Python reserved word)
610
Raymond Hettinger5279fb92011-06-25 11:30:53 +0200611 # we use 1-random() instead of random() to preclude the
612 # possibility of taking the log of zero.
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700613 return -_log(1.0 - self.random()) / lambd
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000614
Tim Petersd7b5e882001-01-25 03:36:26 +0000615 def vonmisesvariate(self, mu, kappa):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000616 """Circular data distribution.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000617
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000618 mu is the mean angle, expressed in radians between 0 and 2*pi, and
619 kappa is the concentration parameter, which must be greater than or
620 equal to zero. If kappa is equal to zero, this distribution reduces
621 to a uniform random angle over the range 0 to 2*pi.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000622
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000623 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000624 # Based upon an algorithm published in: Fisher, N.I.,
625 # "Statistical Analysis of Circular Data", Cambridge
626 # University Press, 1993.
627
628 # Thanks to Magnus Kessler for a correction to the
629 # implementation of step 4.
630
631 random = self.random
632 if kappa <= 1e-6:
633 return TWOPI * random()
634
Serhiy Storchaka6c22b1d2013-02-10 19:28:56 +0200635 s = 0.5 / kappa
636 r = s + _sqrt(1.0 + s * s)
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000637
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700638 while True:
Tim Peters0c9886d2001-01-15 01:18:21 +0000639 u1 = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000640 z = _cos(_pi * u1)
Tim Petersd7b5e882001-01-25 03:36:26 +0000641
Serhiy Storchaka6c22b1d2013-02-10 19:28:56 +0200642 d = z / (r + z)
Tim Petersd7b5e882001-01-25 03:36:26 +0000643 u2 = random()
Serhiy Storchaka6c22b1d2013-02-10 19:28:56 +0200644 if u2 < 1.0 - d * d or u2 <= (1.0 - d) * _exp(d):
Tim Peters0c9886d2001-01-15 01:18:21 +0000645 break
Tim Petersd7b5e882001-01-25 03:36:26 +0000646
Serhiy Storchaka6c22b1d2013-02-10 19:28:56 +0200647 q = 1.0 / r
648 f = (q + z) / (1.0 + q * z)
Tim Petersd7b5e882001-01-25 03:36:26 +0000649 u3 = random()
650 if u3 > 0.5:
Mark Dickinsonbe5f9192013-02-10 14:16:10 +0000651 theta = (mu + _acos(f)) % TWOPI
Tim Petersd7b5e882001-01-25 03:36:26 +0000652 else:
Mark Dickinsonbe5f9192013-02-10 14:16:10 +0000653 theta = (mu - _acos(f)) % TWOPI
Tim Petersd7b5e882001-01-25 03:36:26 +0000654
655 return theta
656
Tim Petersd7b5e882001-01-25 03:36:26 +0000657 def gammavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000658 """Gamma distribution. Not the gamma function!
659
660 Conditions on the parameters are alpha > 0 and beta > 0.
661
Raymond Hettingera8e4d6e2011-03-22 15:55:51 -0700662 The probability distribution function is:
663
664 x ** (alpha - 1) * math.exp(-x / beta)
665 pdf(x) = --------------------------------------
666 math.gamma(alpha) * beta ** alpha
667
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000668 """
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000669 # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
Tim Peters8ac14952002-05-23 15:15:30 +0000670
Guido van Rossum570764d2002-05-14 14:08:12 +0000671 # Warning: a few older sources define the gamma distribution in terms
672 # of alpha > -1.0
673 if alpha <= 0.0 or beta <= 0.0:
Collin Winterce36ad82007-08-30 01:19:48 +0000674 raise ValueError('gammavariate: alpha and beta must be > 0.0')
Tim Peters8ac14952002-05-23 15:15:30 +0000675
Tim Petersd7b5e882001-01-25 03:36:26 +0000676 random = self.random
Tim Petersd7b5e882001-01-25 03:36:26 +0000677 if alpha > 1.0:
678
679 # Uses R.C.H. Cheng, "The generation of Gamma
680 # variables with non-integral shape parameters",
681 # Applied Statistics, (1977), 26, No. 1, p71-74
682
Raymond Hettingerca6cdc22002-05-13 23:40:14 +0000683 ainv = _sqrt(2.0 * alpha - 1.0)
684 bbb = alpha - LOG4
685 ccc = alpha + ainv
Tim Peters8ac14952002-05-23 15:15:30 +0000686
Raymond Hettinger6a613f92020-08-02 12:03:32 -0700687 while True:
Tim Petersd7b5e882001-01-25 03:36:26 +0000688 u1 = random()
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700689 if not 1e-7 < u1 < 0.9999999:
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000690 continue
691 u2 = 1.0 - random()
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700692 v = _log(u1 / (1.0 - u1)) / ainv
693 x = alpha * _exp(v)
694 z = u1 * u1 * u2
695 r = bbb + ccc * v - x
696 if r + SG_MAGICCONST - 4.5 * z >= 0.0 or r >= _log(z):
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000697 return x * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000698
699 elif alpha == 1.0:
leodema9f396b62017-06-04 07:41:41 +0100700 # expovariate(1/beta)
leodema63d15222018-12-24 07:54:25 +0100701 return -_log(1.0 - random()) * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000702
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700703 else:
704 # alpha is between 0 and 1 (exclusive)
Tim Petersd7b5e882001-01-25 03:36:26 +0000705 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700706 while True:
Tim Petersd7b5e882001-01-25 03:36:26 +0000707 u = random()
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700708 b = (_e + alpha) / _e
709 p = b * u
Tim Petersd7b5e882001-01-25 03:36:26 +0000710 if p <= 1.0:
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700711 x = p ** (1.0 / alpha)
Tim Petersd7b5e882001-01-25 03:36:26 +0000712 else:
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700713 x = -_log((b - p) / alpha)
Tim Petersd7b5e882001-01-25 03:36:26 +0000714 u1 = random()
Raymond Hettinger42406e62005-04-30 09:02:51 +0000715 if p > 1.0:
716 if u1 <= x ** (alpha - 1.0):
717 break
718 elif u1 <= _exp(-x):
Tim Petersd7b5e882001-01-25 03:36:26 +0000719 break
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000720 return x * beta
721
Tim Petersd7b5e882001-01-25 03:36:26 +0000722 def betavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000723 """Beta distribution.
724
Thomas Woutersb2137042007-02-01 18:02:27 +0000725 Conditions on the parameters are alpha > 0 and beta > 0.
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000726 Returned values range between 0 and 1.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000727
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000728 """
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700729 ## See
730 ## http://mail.python.org/pipermail/python-bugs-list/2001-January/003752.html
731 ## for Ivan Frohne's insightful analysis of why the original implementation:
732 ##
733 ## def betavariate(self, alpha, beta):
734 ## # Discrete Event Simulation in C, pp 87-88.
735 ##
736 ## y = self.expovariate(alpha)
737 ## z = self.expovariate(1.0/beta)
738 ## return z/(y+z)
739 ##
740 ## was dead wrong, and how it probably got that way.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000741
Tim Peters85e2e472001-01-26 06:49:56 +0000742 # This version due to Janne Sinkkonen, and matches all the std
743 # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
Raymond Hettinger650c1c92016-06-25 05:36:42 +0300744 y = self.gammavariate(alpha, 1.0)
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700745 if y:
Raymond Hettinger650c1c92016-06-25 05:36:42 +0300746 return y / (y + self.gammavariate(beta, 1.0))
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700747 return 0.0
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000748
Tim Petersd7b5e882001-01-25 03:36:26 +0000749 def paretovariate(self, alpha):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000750 """Pareto distribution. alpha is the shape parameter."""
Tim Petersd7b5e882001-01-25 03:36:26 +0000751 # Jain, pg. 495
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000752
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000753 u = 1.0 - self.random()
Raymond Hettinger5c327092020-08-01 01:18:26 -0700754 return u ** (-1.0 / alpha)
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000755
Tim Petersd7b5e882001-01-25 03:36:26 +0000756 def weibullvariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000757 """Weibull distribution.
758
759 alpha is the scale parameter and beta is the shape parameter.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000760
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000761 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000762 # Jain, pg. 499; bug fix courtesy Bill Arms
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000763
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000764 u = 1.0 - self.random()
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700765 return alpha * (-_log(u)) ** (1.0 / beta)
766
Guido van Rossum6c395ba1999-08-18 13:53:28 +0000767
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700768## ------------------------------------------------------------------
Raymond Hettinger23f12412004-09-13 22:23:21 +0000769## --------------- Operating System Random Source ------------------
Raymond Hettinger356a4592004-08-30 06:14:31 +0000770
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700771
Raymond Hettinger23f12412004-09-13 22:23:21 +0000772class SystemRandom(Random):
773 """Alternate random number generator using sources provided
774 by the operating system (such as /dev/urandom on Unix or
775 CryptGenRandom on Windows).
Raymond Hettinger356a4592004-08-30 06:14:31 +0000776
777 Not available on all systems (see os.urandom() for details).
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700778
Raymond Hettinger356a4592004-08-30 06:14:31 +0000779 """
780
781 def random(self):
782 """Get the next random number in the range [0.0, 1.0)."""
Raymond Hettinger183cd1f2010-09-08 18:48:21 +0000783 return (int.from_bytes(_urandom(7), 'big') >> 3) * RECIP_BPF
Raymond Hettinger356a4592004-08-30 06:14:31 +0000784
785 def getrandbits(self, k):
Serhiy Storchaka95949422013-08-27 19:40:23 +0300786 """getrandbits(k) -> x. Generates an int with k random bits."""
Antoine Pitrou75a33782020-04-17 19:32:14 +0200787 if k < 0:
788 raise ValueError('number of bits must be non-negative')
Raymond Hettinger63b17672010-09-08 19:27:59 +0000789 numbytes = (k + 7) // 8 # bits / 8 and rounded up
790 x = int.from_bytes(_urandom(numbytes), 'big')
791 return x >> (numbytes * 8 - k) # trim excess bits
Raymond Hettinger356a4592004-08-30 06:14:31 +0000792
Victor Stinner9f5fe792020-04-17 19:05:35 +0200793 def randbytes(self, n):
794 """Generate n random bytes."""
795 # os.urandom(n) fails with ValueError for n < 0
796 # and returns an empty bytes string for n == 0.
797 return _urandom(n)
798
Raymond Hettinger28de64f2008-01-13 23:40:30 +0000799 def seed(self, *args, **kwds):
Raymond Hettinger23f12412004-09-13 22:23:21 +0000800 "Stub method. Not used for a system random number generator."
Raymond Hettinger356a4592004-08-30 06:14:31 +0000801 return None
Raymond Hettinger356a4592004-08-30 06:14:31 +0000802
803 def _notimplemented(self, *args, **kwds):
Raymond Hettinger23f12412004-09-13 22:23:21 +0000804 "Method should not be called for a system random number generator."
805 raise NotImplementedError('System entropy source does not have state.')
Raymond Hettinger356a4592004-08-30 06:14:31 +0000806 getstate = setstate = _notimplemented
807
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700808
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700809# ----------------------------------------------------------------------
810# Create one instance, seeded from current time, and export its methods
811# as module-level functions. The functions share state across all uses
812# (both in the user's code and in the Python libraries), but that's fine
813# for most programs and is easier for the casual user than making them
814# instantiate their own Random() instance.
815
816_inst = Random()
817seed = _inst.seed
818random = _inst.random
819uniform = _inst.uniform
820triangular = _inst.triangular
821randint = _inst.randint
822choice = _inst.choice
823randrange = _inst.randrange
824sample = _inst.sample
825shuffle = _inst.shuffle
826choices = _inst.choices
827normalvariate = _inst.normalvariate
828lognormvariate = _inst.lognormvariate
829expovariate = _inst.expovariate
830vonmisesvariate = _inst.vonmisesvariate
831gammavariate = _inst.gammavariate
832gauss = _inst.gauss
833betavariate = _inst.betavariate
834paretovariate = _inst.paretovariate
835weibullvariate = _inst.weibullvariate
836getstate = _inst.getstate
837setstate = _inst.setstate
838getrandbits = _inst.getrandbits
839randbytes = _inst.randbytes
840
841
842## ------------------------------------------------------
843## ----------------- test program -----------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000844
Raymond Hettinger62297132003-08-30 01:24:19 +0000845def _test_generator(n, func, args):
Raymond Hettinger26a1ad12020-06-22 19:38:59 -0700846 from statistics import stdev, fmean as mean
847 from time import perf_counter
848
849 t0 = perf_counter()
850 data = [func(*args) for i in range(n)]
851 t1 = perf_counter()
852
853 xbar = mean(data)
854 sigma = stdev(data, xbar)
855 low = min(data)
856 high = max(data)
857
858 print(f'{t1 - t0:.3f} sec, {n} times {func.__name__}')
859 print('avg %g, stddev %g, min %g, max %g\n' % (xbar, sigma, low, high))
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000860
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000861
862def _test(N=2000):
Raymond Hettinger62297132003-08-30 01:24:19 +0000863 _test_generator(N, random, ())
864 _test_generator(N, normalvariate, (0.0, 1.0))
865 _test_generator(N, lognormvariate, (0.0, 1.0))
866 _test_generator(N, vonmisesvariate, (0.0, 1.0))
867 _test_generator(N, gammavariate, (0.01, 1.0))
868 _test_generator(N, gammavariate, (0.1, 1.0))
869 _test_generator(N, gammavariate, (0.1, 2.0))
870 _test_generator(N, gammavariate, (0.5, 1.0))
871 _test_generator(N, gammavariate, (0.9, 1.0))
872 _test_generator(N, gammavariate, (1.0, 1.0))
873 _test_generator(N, gammavariate, (2.0, 1.0))
874 _test_generator(N, gammavariate, (20.0, 1.0))
875 _test_generator(N, gammavariate, (200.0, 1.0))
876 _test_generator(N, gauss, (0.0, 1.0))
877 _test_generator(N, betavariate, (3.0, 3.0))
Raymond Hettinger9db5b8d2020-06-13 09:46:47 -0700878 _test_generator(N, triangular, (0.0, 1.0, 1.0 / 3.0))
Tim Peterscd804102001-01-25 20:25:57 +0000879
Raymond Hettinger40f62172002-12-29 23:03:38 +0000880
Raymond Hettingeref19bad2020-06-25 17:03:50 -0700881## ------------------------------------------------------
882## ------------------ fork support ---------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000883
Antoine Pitrou346cbd32017-05-27 17:50:54 +0200884if hasattr(_os, "fork"):
Gregory P. Smith163468a2017-05-29 10:03:41 -0700885 _os.register_at_fork(after_in_child=_inst.seed)
Antoine Pitrou346cbd32017-05-27 17:50:54 +0200886
887
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000888if __name__ == '__main__':
Tim Petersd7b5e882001-01-25 03:36:26 +0000889 _test()