blob: f2c4f39fb6079d55edf1e331f9adb6a136c5b1c1 [file] [log] [blame]
Guido van Rossume7b146f2000-02-04 15:28:42 +00001"""Random variable generators.
Guido van Rossumff03b1a1994-03-09 12:55:02 +00002
Tim Petersd7b5e882001-01-25 03:36:26 +00003 integers
4 --------
5 uniform within range
6
7 sequences
8 ---------
9 pick random element
Raymond Hettingerf24eb352002-11-12 17:41:57 +000010 pick random sample
Raymond Hettingere8f1e002016-09-06 17:15:29 -070011 pick weighted random sample
Tim Petersd7b5e882001-01-25 03:36:26 +000012 generate random permutation
13
Guido van Rossume7b146f2000-02-04 15:28:42 +000014 distributions on the real line:
15 ------------------------------
Tim Petersd7b5e882001-01-25 03:36:26 +000016 uniform
Christian Heimesfe337bf2008-03-23 21:54:12 +000017 triangular
Guido van Rossume7b146f2000-02-04 15:28:42 +000018 normal (Gaussian)
19 lognormal
20 negative exponential
21 gamma
22 beta
Raymond Hettinger40f62172002-12-29 23:03:38 +000023 pareto
24 Weibull
Guido van Rossumff03b1a1994-03-09 12:55:02 +000025
Guido van Rossume7b146f2000-02-04 15:28:42 +000026 distributions on the circle (angles 0 to 2pi)
27 ---------------------------------------------
28 circular uniform
29 von Mises
30
Raymond Hettinger40f62172002-12-29 23:03:38 +000031General notes on the underlying Mersenne Twister core generator:
Guido van Rossume7b146f2000-02-04 15:28:42 +000032
Raymond Hettinger40f62172002-12-29 23:03:38 +000033* The period is 2**19937-1.
Thomas Wouters0e3f5912006-08-11 14:57:12 +000034* It is one of the most extensively tested generators in existence.
Thomas Wouters0e3f5912006-08-11 14:57:12 +000035* The random() method is implemented in C, executes in a single Python step,
36 and is, therefore, threadsafe.
Tim Peterse360d952001-01-26 10:00:39 +000037
Guido van Rossume7b146f2000-02-04 15:28:42 +000038"""
Guido van Rossumd03e1191998-05-29 17:51:31 +000039
Raymond Hettinger2f726e92003-10-05 09:09:15 +000040from warnings import warn as _warn
Raymond Hettinger91e27c22005-08-19 01:36:35 +000041from math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil
Tim Petersd7b5e882001-01-25 03:36:26 +000042from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
Raymond Hettingerc1c43ca2004-09-05 00:00:42 +000043from os import urandom as _urandom
Christian Heimesf1dc3ee2013-10-13 02:04:20 +020044from _collections_abc import Set as _Set, Sequence as _Sequence
Raymond Hettinger53822032019-02-16 13:30:51 -080045from itertools import accumulate as _accumulate, repeat as _repeat
Raymond Hettingercfd31f02019-02-13 02:04:17 -080046from bisect import bisect as _bisect
Antoine Pitrou346cbd32017-05-27 17:50:54 +020047import os as _os
Guido van Rossumff03b1a1994-03-09 12:55:02 +000048
Christian Heimesd9145962019-04-10 22:18:02 +020049try:
50 # hashlib is pretty heavy to load, try lean internal module first
51 from _sha512 import sha512 as _sha512
52except ImportError:
53 # fallback to official implementation
54 from hashlib import sha512 as _sha512
55
56
Raymond Hettingerf24eb352002-11-12 17:41:57 +000057__all__ = ["Random","seed","random","uniform","randint","choice","sample",
Skip Montanaro0de65802001-02-15 22:15:14 +000058 "randrange","shuffle","normalvariate","lognormvariate",
Christian Heimesfe337bf2008-03-23 21:54:12 +000059 "expovariate","vonmisesvariate","gammavariate","triangular",
Raymond Hettingerf8a52d32003-08-05 12:23:19 +000060 "gauss","betavariate","paretovariate","weibullvariate",
Raymond Hettinger28aa4a02016-09-07 00:08:44 -070061 "getstate","setstate", "getrandbits", "choices",
Raymond Hettinger23f12412004-09-13 22:23:21 +000062 "SystemRandom"]
Tim Petersd7b5e882001-01-25 03:36:26 +000063
64NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
Tim Petersd7b5e882001-01-25 03:36:26 +000065TWOPI = 2.0*_pi
Tim Petersd7b5e882001-01-25 03:36:26 +000066LOG4 = _log(4.0)
Tim Petersd7b5e882001-01-25 03:36:26 +000067SG_MAGICCONST = 1.0 + _log(4.5)
Raymond Hettinger2f726e92003-10-05 09:09:15 +000068BPF = 53 # Number of bits in a float
Tim Peters7c2a85b2004-08-31 02:19:55 +000069RECIP_BPF = 2**-BPF
Tim Petersd7b5e882001-01-25 03:36:26 +000070
Raymond Hettinger356a4592004-08-30 06:14:31 +000071
Tim Petersd7b5e882001-01-25 03:36:26 +000072# Translated by Guido van Rossum from C source provided by
Raymond Hettinger40f62172002-12-29 23:03:38 +000073# Adrian Baddeley. Adapted by Raymond Hettinger for use with
Raymond Hettinger3fa19d72004-08-31 01:05:15 +000074# the Mersenne Twister and os.urandom() core generators.
Tim Petersd7b5e882001-01-25 03:36:26 +000075
Raymond Hettinger145a4a02003-01-07 10:25:55 +000076import _random
Raymond Hettinger40f62172002-12-29 23:03:38 +000077
Raymond Hettinger145a4a02003-01-07 10:25:55 +000078class Random(_random.Random):
Raymond Hettingerc32f0332002-05-23 19:44:49 +000079 """Random number generator base class used by bound module functions.
80
81 Used to instantiate instances of Random to get generators that don't
Raymond Hettinger28de64f2008-01-13 23:40:30 +000082 share state.
Raymond Hettingerc32f0332002-05-23 19:44:49 +000083
84 Class Random can also be subclassed if you want to use a different basic
85 generator of your own devising: in that case, override the following
Raymond Hettinger28de64f2008-01-13 23:40:30 +000086 methods: random(), seed(), getstate(), and setstate().
Benjamin Petersond18de0e2008-07-31 20:21:46 +000087 Optionally, implement a getrandbits() method so that randrange()
Raymond Hettinger2f726e92003-10-05 09:09:15 +000088 can cover arbitrarily large ranges.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +000089
Raymond Hettingerc32f0332002-05-23 19:44:49 +000090 """
Tim Petersd7b5e882001-01-25 03:36:26 +000091
Christian Heimescbf3b5c2007-12-03 21:02:03 +000092 VERSION = 3 # used by getstate/setstate
Tim Petersd7b5e882001-01-25 03:36:26 +000093
94 def __init__(self, x=None):
95 """Initialize an instance.
96
97 Optional argument x controls seeding, as for Random.seed().
98 """
99
100 self.seed(x)
Raymond Hettinger40f62172002-12-29 23:03:38 +0000101 self.gauss_next = None
Tim Petersd7b5e882001-01-25 03:36:26 +0000102
Serhiy Storchaka2085bd02019-06-01 11:00:15 +0300103 def __init_subclass__(cls, /, **kwargs):
Wolfgang Maierba3a87a2018-04-17 17:16:17 +0200104 """Control how subclasses generate random integers.
105
106 The algorithm a subclass can use depends on the random() and/or
107 getrandbits() implementation available to it and determines
108 whether it can generate random integers from arbitrarily large
109 ranges.
110 """
111
Serhiy Storchakaec1622d2018-05-08 15:45:15 +0300112 for c in cls.__mro__:
113 if '_randbelow' in c.__dict__:
114 # just inherit it
115 break
116 if 'getrandbits' in c.__dict__:
117 cls._randbelow = cls._randbelow_with_getrandbits
118 break
119 if 'random' in c.__dict__:
120 cls._randbelow = cls._randbelow_without_getrandbits
121 break
Wolfgang Maierba3a87a2018-04-17 17:16:17 +0200122
Raymond Hettingerf763a722010-09-07 00:38:15 +0000123 def seed(self, a=None, version=2):
Raymond Hettingerd0cdeaa2019-08-22 09:19:36 -0700124 """Initialize internal state from a seed.
125
126 The only supported seed types are None, int, float,
127 str, bytes, and bytearray.
Tim Petersd7b5e882001-01-25 03:36:26 +0000128
Raymond Hettinger23f12412004-09-13 22:23:21 +0000129 None or no argument seeds from current time or from an operating
130 system specific randomness source if available.
Tim Peters0de88fc2001-02-01 04:59:18 +0000131
Raymond Hettinger183cd1f2010-09-08 18:48:21 +0000132 If *a* is an int, all bits are used.
Raymond Hettingerf763a722010-09-07 00:38:15 +0000133
Raymond Hettinger16eb8272016-09-04 11:17:28 -0700134 For version 2 (the default), all of the bits are used if *a* is a str,
135 bytes, or bytearray. For version 1 (provided for reproducing random
136 sequences from older versions of Python), the algorithm for str and
137 bytes generates a narrower range of seeds.
138
Tim Petersd7b5e882001-01-25 03:36:26 +0000139 """
140
Raymond Hettingerc7bab7c2016-08-31 15:01:08 -0700141 if version == 1 and isinstance(a, (str, bytes)):
Raymond Hettinger132a7d72017-09-17 09:04:30 -0700142 a = a.decode('latin-1') if isinstance(a, bytes) else a
Raymond Hettingerc7bab7c2016-08-31 15:01:08 -0700143 x = ord(a[0]) << 7 if a else 0
Raymond Hettinger132a7d72017-09-17 09:04:30 -0700144 for c in map(ord, a):
145 x = ((1000003 * x) ^ c) & 0xFFFFFFFFFFFFFFFF
Raymond Hettingerc7bab7c2016-08-31 15:01:08 -0700146 x ^= len(a)
147 a = -2 if x == -1 else x
148
Raymond Hettingerd0cdeaa2019-08-22 09:19:36 -0700149 elif version == 2 and isinstance(a, (str, bytes, bytearray)):
Raymond Hettinger2f9cc7a2016-08-31 23:00:32 -0700150 if isinstance(a, str):
151 a = a.encode()
152 a += _sha512(a).digest()
153 a = int.from_bytes(a, 'big')
Raymond Hettingerf763a722010-09-07 00:38:15 +0000154
Raymond Hettingerd0cdeaa2019-08-22 09:19:36 -0700155 elif not isinstance(a, (type(None), int, float, str, bytes, bytearray)):
156 _warn('Seeding based on hashing is deprecated\n'
157 'since Python 3.9 and will be removed in a subsequent '
158 'version. The only \n'
159 'supported seed types are: None, '
160 'int, float, str, bytes, and bytearray.',
161 DeprecationWarning, 2)
162
Guido van Rossumcd16bf62007-06-13 18:07:49 +0000163 super().seed(a)
Tim Peters46c04e12002-05-05 20:40:00 +0000164 self.gauss_next = None
165
Tim Peterscd804102001-01-25 20:25:57 +0000166 def getstate(self):
167 """Return internal state; can be passed to setstate() later."""
Guido van Rossumcd16bf62007-06-13 18:07:49 +0000168 return self.VERSION, super().getstate(), self.gauss_next
Tim Peterscd804102001-01-25 20:25:57 +0000169
170 def setstate(self, state):
171 """Restore internal state from object returned by getstate()."""
172 version = state[0]
Christian Heimescbf3b5c2007-12-03 21:02:03 +0000173 if version == 3:
Raymond Hettinger40f62172002-12-29 23:03:38 +0000174 version, internalstate, self.gauss_next = state
Guido van Rossumcd16bf62007-06-13 18:07:49 +0000175 super().setstate(internalstate)
Christian Heimescbf3b5c2007-12-03 21:02:03 +0000176 elif version == 2:
177 version, internalstate, self.gauss_next = state
178 # In version 2, the state was saved as signed ints, which causes
179 # inconsistencies between 32/64-bit systems. The state is
180 # really unsigned 32-bit ints, so we convert negative ints from
181 # version 2 to positive longs for version 3.
182 try:
Raymond Hettingerc585eec2010-09-07 15:00:15 +0000183 internalstate = tuple(x % (2**32) for x in internalstate)
Christian Heimescbf3b5c2007-12-03 21:02:03 +0000184 except ValueError as e:
185 raise TypeError from e
Raymond Hettinger183cd1f2010-09-08 18:48:21 +0000186 super().setstate(internalstate)
Tim Peterscd804102001-01-25 20:25:57 +0000187 else:
188 raise ValueError("state with version %s passed to "
189 "Random.setstate() of version %s" %
190 (version, self.VERSION))
191
Tim Peterscd804102001-01-25 20:25:57 +0000192## ---- Methods below this point do not need to be overridden when
193## ---- subclassing for the purpose of using a different core generator.
194
Victor Stinner2d875772020-04-29 18:49:00 +0200195## -------------------- bytes methods ---------------------
196
197 def randbytes(self, n):
198 """Generate n random bytes."""
199 return self.getrandbits(n * 8).to_bytes(n, 'little')
200
Tim Peterscd804102001-01-25 20:25:57 +0000201## -------------------- pickle support -------------------
202
R David Murrayd9ebf4d2013-04-02 13:10:52 -0400203 # Issue 17489: Since __reduce__ was defined to fix #759889 this is no
204 # longer called; we leave it here because it has been here since random was
205 # rewritten back in 2001 and why risk breaking something.
Tim Peterscd804102001-01-25 20:25:57 +0000206 def __getstate__(self): # for pickle
207 return self.getstate()
208
209 def __setstate__(self, state): # for pickle
210 self.setstate(state)
211
Raymond Hettinger5f078ff2003-06-24 20:29:04 +0000212 def __reduce__(self):
213 return self.__class__, (), self.getstate()
214
Tim Peterscd804102001-01-25 20:25:57 +0000215## -------------------- integer methods -------------------
216
Raymond Hettinger8fe47c32013-10-05 21:48:21 -0700217 def randrange(self, start, stop=None, step=1, _int=int):
Tim Petersd7b5e882001-01-25 03:36:26 +0000218 """Choose a random item from range(start, stop[, step]).
219
220 This fixes the problem with randint() which includes the
221 endpoint; in Python this is usually not what you want.
Raymond Hettinger3051cc32010-09-07 00:48:40 +0000222
Tim Petersd7b5e882001-01-25 03:36:26 +0000223 """
224
225 # This code is a bit messy to make it fast for the
Tim Peters9146f272002-08-16 03:41:39 +0000226 # common case while still doing adequate error checking.
Raymond Hettinger8fe47c32013-10-05 21:48:21 -0700227 istart = _int(start)
Tim Petersd7b5e882001-01-25 03:36:26 +0000228 if istart != start:
Collin Winterce36ad82007-08-30 01:19:48 +0000229 raise ValueError("non-integer arg 1 for randrange()")
Raymond Hettinger3051cc32010-09-07 00:48:40 +0000230 if stop is None:
Tim Petersd7b5e882001-01-25 03:36:26 +0000231 if istart > 0:
Raymond Hettinger05156612010-09-07 04:44:52 +0000232 return self._randbelow(istart)
Collin Winterce36ad82007-08-30 01:19:48 +0000233 raise ValueError("empty range for randrange()")
Tim Peters9146f272002-08-16 03:41:39 +0000234
235 # stop argument supplied.
Raymond Hettinger8fe47c32013-10-05 21:48:21 -0700236 istop = _int(stop)
Tim Petersd7b5e882001-01-25 03:36:26 +0000237 if istop != stop:
Collin Winterce36ad82007-08-30 01:19:48 +0000238 raise ValueError("non-integer stop for randrange()")
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000239 width = istop - istart
240 if step == 1 and width > 0:
Raymond Hettingerc3246972010-09-07 09:32:57 +0000241 return istart + self._randbelow(width)
Tim Petersd7b5e882001-01-25 03:36:26 +0000242 if step == 1:
Kumar Akshay2433a2a2019-01-22 00:49:59 +0530243 raise ValueError("empty range for randrange() (%d, %d, %d)" % (istart, istop, width))
Tim Peters9146f272002-08-16 03:41:39 +0000244
245 # Non-unit step argument supplied.
Raymond Hettinger8fe47c32013-10-05 21:48:21 -0700246 istep = _int(step)
Tim Petersd7b5e882001-01-25 03:36:26 +0000247 if istep != step:
Collin Winterce36ad82007-08-30 01:19:48 +0000248 raise ValueError("non-integer step for randrange()")
Tim Petersd7b5e882001-01-25 03:36:26 +0000249 if istep > 0:
Raymond Hettingerffdb8bb2004-09-27 15:29:05 +0000250 n = (width + istep - 1) // istep
Tim Petersd7b5e882001-01-25 03:36:26 +0000251 elif istep < 0:
Raymond Hettingerffdb8bb2004-09-27 15:29:05 +0000252 n = (width + istep + 1) // istep
Tim Petersd7b5e882001-01-25 03:36:26 +0000253 else:
Collin Winterce36ad82007-08-30 01:19:48 +0000254 raise ValueError("zero step for randrange()")
Tim Petersd7b5e882001-01-25 03:36:26 +0000255
256 if n <= 0:
Collin Winterce36ad82007-08-30 01:19:48 +0000257 raise ValueError("empty range for randrange()")
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000258
Raymond Hettinger05156612010-09-07 04:44:52 +0000259 return istart + istep*self._randbelow(n)
Tim Petersd7b5e882001-01-25 03:36:26 +0000260
261 def randint(self, a, b):
Tim Peterscd804102001-01-25 20:25:57 +0000262 """Return random integer in range [a, b], including both end points.
Tim Petersd7b5e882001-01-25 03:36:26 +0000263 """
264
265 return self.randrange(a, b+1)
266
Wolfgang Maierba3a87a2018-04-17 17:16:17 +0200267 def _randbelow_with_getrandbits(self, n):
Raymond Hettinger4168f1e2020-05-01 10:34:19 -0700268 "Return a random int in the range [0,n). Returns 0 if n==0."
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000269
Antoine Pitrou75a33782020-04-17 19:32:14 +0200270 if not n:
Raymond Hettinger4168f1e2020-05-01 10:34:19 -0700271 return 0
Raymond Hettingerc3246972010-09-07 09:32:57 +0000272 getrandbits = self.getrandbits
Wolfgang Maierba3a87a2018-04-17 17:16:17 +0200273 k = n.bit_length() # don't use (n-1) here because n can be 1
274 r = getrandbits(k) # 0 <= r < 2**k
275 while r >= n:
276 r = getrandbits(k)
277 return r
278
279 def _randbelow_without_getrandbits(self, n, int=int, maxsize=1<<BPF):
Raymond Hettinger4168f1e2020-05-01 10:34:19 -0700280 """Return a random int in the range [0,n). Returns 0 if n==0.
Wolfgang Maierba3a87a2018-04-17 17:16:17 +0200281
282 The implementation does not use getrandbits, but only random.
283 """
284
285 random = self.random
Raymond Hettingere4a3e992010-09-08 00:30:28 +0000286 if n >= maxsize:
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000287 _warn("Underlying random() generator does not supply \n"
Raymond Hettingerf015b3f2010-09-07 20:04:42 +0000288 "enough bits to choose from a population range this large.\n"
289 "To remove the range limitation, add a getrandbits() method.")
Raymond Hettingere4a3e992010-09-08 00:30:28 +0000290 return int(random() * n)
Wolfgang Maier091e95e2018-04-05 17:19:44 +0200291 if n == 0:
Raymond Hettinger4168f1e2020-05-01 10:34:19 -0700292 return 0
Raymond Hettingere4a3e992010-09-08 00:30:28 +0000293 rem = maxsize % n
294 limit = (maxsize - rem) / maxsize # int(limit * maxsize) % n == 0
295 r = random()
296 while r >= limit:
297 r = random()
298 return int(r*maxsize) % n
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000299
Wolfgang Maierba3a87a2018-04-17 17:16:17 +0200300 _randbelow = _randbelow_with_getrandbits
301
Tim Peterscd804102001-01-25 20:25:57 +0000302## -------------------- sequence methods -------------------
303
Tim Petersd7b5e882001-01-25 03:36:26 +0000304 def choice(self, seq):
305 """Choose a random element from a non-empty sequence."""
Raymond Hettinger4168f1e2020-05-01 10:34:19 -0700306 return seq[self._randbelow(len(seq))] # raises IndexError if seq is empty
Tim Petersd7b5e882001-01-25 03:36:26 +0000307
Raymond Hettinger8fe47c32013-10-05 21:48:21 -0700308 def shuffle(self, x, random=None):
Antoine Pitrou5e394332012-11-04 02:10:33 +0100309 """Shuffle list x in place, and return None.
Tim Petersd7b5e882001-01-25 03:36:26 +0000310
Antoine Pitrou5e394332012-11-04 02:10:33 +0100311 Optional argument random is a 0-argument function returning a
312 random float in [0.0, 1.0); if it is the default None, the
313 standard random.random will be used.
Senthil Kumaranf8ce51a2013-09-11 22:54:31 -0700314
Tim Petersd7b5e882001-01-25 03:36:26 +0000315 """
316
Raymond Hettinger8fe47c32013-10-05 21:48:21 -0700317 if random is None:
318 randbelow = self._randbelow
319 for i in reversed(range(1, len(x))):
320 # pick an element in x[:i+1] with which to exchange x[i]
321 j = randbelow(i+1)
322 x[i], x[j] = x[j], x[i]
323 else:
Raymond Hettinger190fac92020-05-02 16:45:32 -0700324 _warn('The *random* parameter to shuffle() has been deprecated\n'
325 'since Python 3.9 and will be removed in a subsequent '
326 'version.',
327 DeprecationWarning, 2)
Raymond Hettinger8fe47c32013-10-05 21:48:21 -0700328 _int = int
329 for i in reversed(range(1, len(x))):
330 # pick an element in x[:i+1] with which to exchange x[i]
331 j = _int(random() * (i+1))
332 x[i], x[j] = x[j], x[i]
Tim Petersd7b5e882001-01-25 03:36:26 +0000333
Raymond Hettingerfdbe5222003-06-13 07:01:51 +0000334 def sample(self, population, k):
Raymond Hettinger1acde192008-01-14 01:00:53 +0000335 """Chooses k unique random elements from a population sequence or set.
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000336
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000337 Returns a new list containing elements from the population while
338 leaving the original population unchanged. The resulting list is
339 in selection order so that all sub-slices will also be valid random
340 samples. This allows raffle winners (the sample) to be partitioned
341 into grand prize and second place winners (the subslices).
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000342
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000343 Members of the population need not be hashable or unique. If the
344 population contains repeats, then each occurrence is a possible
345 selection in the sample.
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000346
Guido van Rossum805365e2007-05-07 22:24:25 +0000347 To choose a sample in a range of integers, use range as an argument.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000348 This is especially fast and space efficient for sampling from a
Guido van Rossum805365e2007-05-07 22:24:25 +0000349 large population: sample(range(10000000), 60)
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000350 """
351
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000352 # Sampling without replacement entails tracking either potential
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000353 # selections (the pool) in a list or previous selections in a set.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000354
Jeremy Hylton2b55d352004-02-23 17:27:57 +0000355 # When the number of selections is small compared to the
356 # population, then tracking selections is efficient, requiring
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000357 # only a small set and an occasional reselection. For
Jeremy Hylton2b55d352004-02-23 17:27:57 +0000358 # a larger number of selections, the pool tracking method is
359 # preferred since the list takes less space than the
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000360 # set and it doesn't suffer from frequent reselections.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000361
Raymond Hettinger7fc633f2018-12-04 00:13:38 -0800362 # The number of calls to _randbelow() is kept at or near k, the
363 # theoretical minimum. This is important because running time
364 # is dominated by _randbelow() and because it extracts the
365 # least entropy from the underlying random number generators.
366
367 # Memory requirements are kept to the smaller of a k-length
368 # set or an n-length list.
369
370 # There are other sampling algorithms that do not require
371 # auxiliary memory, but they were rejected because they made
372 # too many calls to _randbelow(), making them slower and
373 # causing them to eat more entropy than necessary.
374
Raymond Hettinger57d1a882011-02-23 00:46:28 +0000375 if isinstance(population, _Set):
Raymond Hettinger4fe00202020-04-19 00:36:42 -0700376 _warn('Sampling from a set deprecated\n'
377 'since Python 3.9 and will be removed in a subsequent version.',
378 DeprecationWarning, 2)
Raymond Hettinger1acde192008-01-14 01:00:53 +0000379 population = tuple(population)
Raymond Hettinger57d1a882011-02-23 00:46:28 +0000380 if not isinstance(population, _Sequence):
Raymond Hettinger4fe00202020-04-19 00:36:42 -0700381 raise TypeError("Population must be a sequence. For dicts or sets, use sorted(d).")
Raymond Hettinger05a505f2010-09-07 19:19:33 +0000382 randbelow = self._randbelow
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000383 n = len(population)
384 if not 0 <= k <= n:
Raymond Hettingerbf871262016-11-21 14:34:33 -0800385 raise ValueError("Sample larger than population or is negative")
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000386 result = [None] * k
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000387 setsize = 21 # size of a small set minus size of an empty list
388 if k > 5:
Tim Peters9e34c042005-08-26 15:20:46 +0000389 setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets
Raymond Hettinger1acde192008-01-14 01:00:53 +0000390 if n <= setsize:
391 # An n-length list is smaller than a k-length set
Raymond Hettinger311f4192002-11-18 09:01:24 +0000392 pool = list(population)
Guido van Rossum805365e2007-05-07 22:24:25 +0000393 for i in range(k): # invariant: non-selected at [0,n-i)
Raymond Hettinger05a505f2010-09-07 19:19:33 +0000394 j = randbelow(n-i)
Raymond Hettinger311f4192002-11-18 09:01:24 +0000395 result[i] = pool[j]
Raymond Hettinger8b9aa8d2003-01-04 05:20:33 +0000396 pool[j] = pool[n-i-1] # move non-selected item into vacancy
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000397 else:
Raymond Hettinger1acde192008-01-14 01:00:53 +0000398 selected = set()
399 selected_add = selected.add
400 for i in range(k):
Raymond Hettinger05a505f2010-09-07 19:19:33 +0000401 j = randbelow(n)
Raymond Hettinger1acde192008-01-14 01:00:53 +0000402 while j in selected:
Raymond Hettinger05a505f2010-09-07 19:19:33 +0000403 j = randbelow(n)
Raymond Hettinger1acde192008-01-14 01:00:53 +0000404 selected_add(j)
405 result[i] = population[j]
Raymond Hettinger311f4192002-11-18 09:01:24 +0000406 return result
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000407
Raymond Hettinger9016f282016-09-26 21:45:57 -0700408 def choices(self, population, weights=None, *, cum_weights=None, k=1):
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700409 """Return a k sized list of population elements chosen with replacement.
410
411 If the relative weights or cumulative weights are not specified,
412 the selections are made with equal probability.
413
414 """
Raymond Hettinger30d00e52016-10-29 16:55:36 -0700415 random = self.random
Raymond Hettingere69cd162018-07-04 15:28:20 -0700416 n = len(population)
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700417 if cum_weights is None:
418 if weights is None:
Raymond Hettinger30d00e52016-10-29 16:55:36 -0700419 _int = int
Raymond Hettinger0a18e052018-11-09 02:39:50 -0800420 n += 0.0 # convert to float for a small speed improvement
Raymond Hettinger53822032019-02-16 13:30:51 -0800421 return [population[_int(random() * n)] for i in _repeat(None, k)]
Raymond Hettingercfd31f02019-02-13 02:04:17 -0800422 cum_weights = list(_accumulate(weights))
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700423 elif weights is not None:
Raymond Hettinger24e42392016-11-13 00:42:56 -0500424 raise TypeError('Cannot specify both weights and cumulative weights')
Raymond Hettingere69cd162018-07-04 15:28:20 -0700425 if len(cum_weights) != n:
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700426 raise ValueError('The number of weights does not match the population')
Raymond Hettinger0a18e052018-11-09 02:39:50 -0800427 total = cum_weights[-1] + 0.0 # convert to float
Raymond Hettinger041d8b42019-11-23 02:22:13 -0800428 if total <= 0.0:
429 raise ValueError('Total of weights must be greater than zero')
430 bisect = _bisect
Raymond Hettingere69cd162018-07-04 15:28:20 -0700431 hi = n - 1
Raymond Hettingerddf71712018-06-27 01:08:31 -0700432 return [population[bisect(cum_weights, random() * total, 0, hi)]
Raymond Hettinger53822032019-02-16 13:30:51 -0800433 for i in _repeat(None, k)]
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700434
Tim Peterscd804102001-01-25 20:25:57 +0000435## -------------------- real-valued distributions -------------------
436
437## -------------------- uniform distribution -------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000438
439 def uniform(self, a, b):
Raymond Hettingerbe40db02009-06-11 23:12:14 +0000440 "Get a random number in the range [a, b) or [a, b] depending on rounding."
Tim Petersd7b5e882001-01-25 03:36:26 +0000441 return a + (b-a) * self.random()
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000442
Christian Heimesfe337bf2008-03-23 21:54:12 +0000443## -------------------- triangular --------------------
444
445 def triangular(self, low=0.0, high=1.0, mode=None):
446 """Triangular distribution.
447
448 Continuous distribution bounded by given lower and upper limits,
449 and having a given mode value in-between.
450
451 http://en.wikipedia.org/wiki/Triangular_distribution
452
453 """
454 u = self.random()
Raymond Hettinger978c6ab2014-05-25 17:25:27 -0700455 try:
456 c = 0.5 if mode is None else (mode - low) / (high - low)
457 except ZeroDivisionError:
458 return low
Christian Heimesfe337bf2008-03-23 21:54:12 +0000459 if u > c:
460 u = 1.0 - u
461 c = 1.0 - c
462 low, high = high, low
Raymond Hettingerf5ea83f2017-09-04 16:51:06 -0700463 return low + (high - low) * _sqrt(u * c)
Christian Heimesfe337bf2008-03-23 21:54:12 +0000464
Tim Peterscd804102001-01-25 20:25:57 +0000465## -------------------- normal distribution --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000466
Tim Petersd7b5e882001-01-25 03:36:26 +0000467 def normalvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000468 """Normal distribution.
469
470 mu is the mean, and sigma is the standard deviation.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000471
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000472 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000473 # mu = mean, sigma = standard deviation
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000474
Tim Petersd7b5e882001-01-25 03:36:26 +0000475 # Uses Kinderman and Monahan method. Reference: Kinderman,
476 # A.J. and Monahan, J.F., "Computer generation of random
477 # variables using the ratio of uniform deviates", ACM Trans
478 # Math Software, 3, (1977), pp257-260.
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000479
Tim Petersd7b5e882001-01-25 03:36:26 +0000480 random = self.random
Raymond Hettinger42406e62005-04-30 09:02:51 +0000481 while 1:
Tim Peters0c9886d2001-01-15 01:18:21 +0000482 u1 = random()
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000483 u2 = 1.0 - random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000484 z = NV_MAGICCONST*(u1-0.5)/u2
485 zz = z*z/4.0
486 if zz <= -_log(u2):
487 break
488 return mu + z*sigma
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000489
Tim Peterscd804102001-01-25 20:25:57 +0000490## -------------------- lognormal distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000491
492 def lognormvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000493 """Log normal distribution.
494
495 If you take the natural logarithm of this distribution, you'll get a
496 normal distribution with mean mu and standard deviation sigma.
497 mu can have any value, and sigma must be greater than zero.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000498
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000499 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000500 return _exp(self.normalvariate(mu, sigma))
501
Tim Peterscd804102001-01-25 20:25:57 +0000502## -------------------- exponential distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000503
504 def expovariate(self, lambd):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000505 """Exponential distribution.
506
Mark Dickinson2f947362009-01-07 17:54:07 +0000507 lambd is 1.0 divided by the desired mean. It should be
508 nonzero. (The parameter would be called "lambda", but that is
509 a reserved word in Python.) Returned values range from 0 to
510 positive infinity if lambd is positive, and from negative
511 infinity to 0 if lambd is negative.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000512
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000513 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000514 # lambd: rate lambd = 1/mean
515 # ('lambda' is a Python reserved word)
516
Raymond Hettinger5279fb92011-06-25 11:30:53 +0200517 # we use 1-random() instead of random() to preclude the
518 # possibility of taking the log of zero.
519 return -_log(1.0 - self.random())/lambd
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000520
Tim Peterscd804102001-01-25 20:25:57 +0000521## -------------------- von Mises distribution --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000522
Tim Petersd7b5e882001-01-25 03:36:26 +0000523 def vonmisesvariate(self, mu, kappa):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000524 """Circular data distribution.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000525
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000526 mu is the mean angle, expressed in radians between 0 and 2*pi, and
527 kappa is the concentration parameter, which must be greater than or
528 equal to zero. If kappa is equal to zero, this distribution reduces
529 to a uniform random angle over the range 0 to 2*pi.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000530
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000531 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000532 # mu: mean angle (in radians between 0 and 2*pi)
533 # kappa: concentration parameter kappa (>= 0)
534 # if kappa = 0 generate uniform random angle
535
536 # Based upon an algorithm published in: Fisher, N.I.,
537 # "Statistical Analysis of Circular Data", Cambridge
538 # University Press, 1993.
539
540 # Thanks to Magnus Kessler for a correction to the
541 # implementation of step 4.
542
543 random = self.random
544 if kappa <= 1e-6:
545 return TWOPI * random()
546
Serhiy Storchaka6c22b1d2013-02-10 19:28:56 +0200547 s = 0.5 / kappa
548 r = s + _sqrt(1.0 + s * s)
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000549
Raymond Hettinger42406e62005-04-30 09:02:51 +0000550 while 1:
Tim Peters0c9886d2001-01-15 01:18:21 +0000551 u1 = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000552 z = _cos(_pi * u1)
Tim Petersd7b5e882001-01-25 03:36:26 +0000553
Serhiy Storchaka6c22b1d2013-02-10 19:28:56 +0200554 d = z / (r + z)
Tim Petersd7b5e882001-01-25 03:36:26 +0000555 u2 = random()
Serhiy Storchaka6c22b1d2013-02-10 19:28:56 +0200556 if u2 < 1.0 - d * d or u2 <= (1.0 - d) * _exp(d):
Tim Peters0c9886d2001-01-15 01:18:21 +0000557 break
Tim Petersd7b5e882001-01-25 03:36:26 +0000558
Serhiy Storchaka6c22b1d2013-02-10 19:28:56 +0200559 q = 1.0 / r
560 f = (q + z) / (1.0 + q * z)
Tim Petersd7b5e882001-01-25 03:36:26 +0000561 u3 = random()
562 if u3 > 0.5:
Mark Dickinsonbe5f9192013-02-10 14:16:10 +0000563 theta = (mu + _acos(f)) % TWOPI
Tim Petersd7b5e882001-01-25 03:36:26 +0000564 else:
Mark Dickinsonbe5f9192013-02-10 14:16:10 +0000565 theta = (mu - _acos(f)) % TWOPI
Tim Petersd7b5e882001-01-25 03:36:26 +0000566
567 return theta
568
Tim Peterscd804102001-01-25 20:25:57 +0000569## -------------------- gamma distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000570
571 def gammavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000572 """Gamma distribution. Not the gamma function!
573
574 Conditions on the parameters are alpha > 0 and beta > 0.
575
Raymond Hettingera8e4d6e2011-03-22 15:55:51 -0700576 The probability distribution function is:
577
578 x ** (alpha - 1) * math.exp(-x / beta)
579 pdf(x) = --------------------------------------
580 math.gamma(alpha) * beta ** alpha
581
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000582 """
Tim Peters8ac14952002-05-23 15:15:30 +0000583
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000584 # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
Tim Peters8ac14952002-05-23 15:15:30 +0000585
Guido van Rossum570764d2002-05-14 14:08:12 +0000586 # Warning: a few older sources define the gamma distribution in terms
587 # of alpha > -1.0
588 if alpha <= 0.0 or beta <= 0.0:
Collin Winterce36ad82007-08-30 01:19:48 +0000589 raise ValueError('gammavariate: alpha and beta must be > 0.0')
Tim Peters8ac14952002-05-23 15:15:30 +0000590
Tim Petersd7b5e882001-01-25 03:36:26 +0000591 random = self.random
Tim Petersd7b5e882001-01-25 03:36:26 +0000592 if alpha > 1.0:
593
594 # Uses R.C.H. Cheng, "The generation of Gamma
595 # variables with non-integral shape parameters",
596 # Applied Statistics, (1977), 26, No. 1, p71-74
597
Raymond Hettingerca6cdc22002-05-13 23:40:14 +0000598 ainv = _sqrt(2.0 * alpha - 1.0)
599 bbb = alpha - LOG4
600 ccc = alpha + ainv
Tim Peters8ac14952002-05-23 15:15:30 +0000601
Raymond Hettinger42406e62005-04-30 09:02:51 +0000602 while 1:
Tim Petersd7b5e882001-01-25 03:36:26 +0000603 u1 = random()
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000604 if not 1e-7 < u1 < .9999999:
605 continue
606 u2 = 1.0 - random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000607 v = _log(u1/(1.0-u1))/ainv
608 x = alpha*_exp(v)
609 z = u1*u1*u2
610 r = bbb+ccc*v-x
611 if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z):
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000612 return x * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000613
614 elif alpha == 1.0:
leodema9f396b62017-06-04 07:41:41 +0100615 # expovariate(1/beta)
leodema63d15222018-12-24 07:54:25 +0100616 return -_log(1.0 - random()) * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000617
618 else: # alpha is between 0 and 1 (exclusive)
619
620 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
621
Raymond Hettinger42406e62005-04-30 09:02:51 +0000622 while 1:
Tim Petersd7b5e882001-01-25 03:36:26 +0000623 u = random()
624 b = (_e + alpha)/_e
625 p = b*u
626 if p <= 1.0:
Raymond Hettinger42406e62005-04-30 09:02:51 +0000627 x = p ** (1.0/alpha)
Tim Petersd7b5e882001-01-25 03:36:26 +0000628 else:
Tim Petersd7b5e882001-01-25 03:36:26 +0000629 x = -_log((b-p)/alpha)
630 u1 = random()
Raymond Hettinger42406e62005-04-30 09:02:51 +0000631 if p > 1.0:
632 if u1 <= x ** (alpha - 1.0):
633 break
634 elif u1 <= _exp(-x):
Tim Petersd7b5e882001-01-25 03:36:26 +0000635 break
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000636 return x * beta
637
Tim Peterscd804102001-01-25 20:25:57 +0000638## -------------------- Gauss (faster alternative) --------------------
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000639
Tim Petersd7b5e882001-01-25 03:36:26 +0000640 def gauss(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000641 """Gaussian distribution.
642
643 mu is the mean, and sigma is the standard deviation. This is
644 slightly faster than the normalvariate() function.
645
646 Not thread-safe without a lock around calls.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000647
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000648 """
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000649
Tim Petersd7b5e882001-01-25 03:36:26 +0000650 # When x and y are two variables from [0, 1), uniformly
651 # distributed, then
652 #
653 # cos(2*pi*x)*sqrt(-2*log(1-y))
654 # sin(2*pi*x)*sqrt(-2*log(1-y))
655 #
656 # are two *independent* variables with normal distribution
657 # (mu = 0, sigma = 1).
658 # (Lambert Meertens)
659 # (corrected version; bug discovered by Mike Miller, fixed by LM)
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000660
Tim Petersd7b5e882001-01-25 03:36:26 +0000661 # Multithreading note: When two threads call this function
662 # simultaneously, it is possible that they will receive the
663 # same return value. The window is very small though. To
664 # avoid this, you have to use a lock around all calls. (I
665 # didn't want to slow this down in the serial case by using a
666 # lock here.)
Guido van Rossumd03e1191998-05-29 17:51:31 +0000667
Tim Petersd7b5e882001-01-25 03:36:26 +0000668 random = self.random
669 z = self.gauss_next
670 self.gauss_next = None
671 if z is None:
672 x2pi = random() * TWOPI
673 g2rad = _sqrt(-2.0 * _log(1.0 - random()))
674 z = _cos(x2pi) * g2rad
675 self.gauss_next = _sin(x2pi) * g2rad
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000676
Tim Petersd7b5e882001-01-25 03:36:26 +0000677 return mu + z*sigma
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000678
Tim Peterscd804102001-01-25 20:25:57 +0000679## -------------------- beta --------------------
Tim Peters85e2e472001-01-26 06:49:56 +0000680## See
Ezio Melotti20f53f12011-04-15 08:25:16 +0300681## http://mail.python.org/pipermail/python-bugs-list/2001-January/003752.html
Tim Peters85e2e472001-01-26 06:49:56 +0000682## for Ivan Frohne's insightful analysis of why the original implementation:
683##
684## def betavariate(self, alpha, beta):
685## # Discrete Event Simulation in C, pp 87-88.
686##
687## y = self.expovariate(alpha)
688## z = self.expovariate(1.0/beta)
689## return z/(y+z)
690##
691## was dead wrong, and how it probably got that way.
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000692
Tim Petersd7b5e882001-01-25 03:36:26 +0000693 def betavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000694 """Beta distribution.
695
Thomas Woutersb2137042007-02-01 18:02:27 +0000696 Conditions on the parameters are alpha > 0 and beta > 0.
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000697 Returned values range between 0 and 1.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000698
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000699 """
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000700
Tim Peters85e2e472001-01-26 06:49:56 +0000701 # This version due to Janne Sinkkonen, and matches all the std
702 # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
Raymond Hettinger650c1c92016-06-25 05:36:42 +0300703 y = self.gammavariate(alpha, 1.0)
Tim Peters85e2e472001-01-26 06:49:56 +0000704 if y == 0:
705 return 0.0
706 else:
Raymond Hettinger650c1c92016-06-25 05:36:42 +0300707 return y / (y + self.gammavariate(beta, 1.0))
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000708
Tim Peterscd804102001-01-25 20:25:57 +0000709## -------------------- Pareto --------------------
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000710
Tim Petersd7b5e882001-01-25 03:36:26 +0000711 def paretovariate(self, alpha):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000712 """Pareto distribution. alpha is the shape parameter."""
Tim Petersd7b5e882001-01-25 03:36:26 +0000713 # Jain, pg. 495
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000714
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000715 u = 1.0 - self.random()
Raymond Hettinger8ff10992010-09-08 18:58:33 +0000716 return 1.0 / u ** (1.0/alpha)
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000717
Tim Peterscd804102001-01-25 20:25:57 +0000718## -------------------- Weibull --------------------
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000719
Tim Petersd7b5e882001-01-25 03:36:26 +0000720 def weibullvariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000721 """Weibull distribution.
722
723 alpha is the scale parameter and beta is the shape parameter.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000724
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000725 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000726 # Jain, pg. 499; bug fix courtesy Bill Arms
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000727
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000728 u = 1.0 - self.random()
Raymond Hettinger183cd1f2010-09-08 18:48:21 +0000729 return alpha * (-_log(u)) ** (1.0/beta)
Guido van Rossum6c395ba1999-08-18 13:53:28 +0000730
Raymond Hettinger23f12412004-09-13 22:23:21 +0000731## --------------- Operating System Random Source ------------------
Raymond Hettinger356a4592004-08-30 06:14:31 +0000732
Raymond Hettinger23f12412004-09-13 22:23:21 +0000733class SystemRandom(Random):
734 """Alternate random number generator using sources provided
735 by the operating system (such as /dev/urandom on Unix or
736 CryptGenRandom on Windows).
Raymond Hettinger356a4592004-08-30 06:14:31 +0000737
738 Not available on all systems (see os.urandom() for details).
739 """
740
741 def random(self):
742 """Get the next random number in the range [0.0, 1.0)."""
Raymond Hettinger183cd1f2010-09-08 18:48:21 +0000743 return (int.from_bytes(_urandom(7), 'big') >> 3) * RECIP_BPF
Raymond Hettinger356a4592004-08-30 06:14:31 +0000744
745 def getrandbits(self, k):
Serhiy Storchaka95949422013-08-27 19:40:23 +0300746 """getrandbits(k) -> x. Generates an int with k random bits."""
Antoine Pitrou75a33782020-04-17 19:32:14 +0200747 if k < 0:
748 raise ValueError('number of bits must be non-negative')
Raymond Hettinger63b17672010-09-08 19:27:59 +0000749 numbytes = (k + 7) // 8 # bits / 8 and rounded up
750 x = int.from_bytes(_urandom(numbytes), 'big')
751 return x >> (numbytes * 8 - k) # trim excess bits
Raymond Hettinger356a4592004-08-30 06:14:31 +0000752
Victor Stinner9f5fe792020-04-17 19:05:35 +0200753 def randbytes(self, n):
754 """Generate n random bytes."""
755 # os.urandom(n) fails with ValueError for n < 0
756 # and returns an empty bytes string for n == 0.
757 return _urandom(n)
758
Raymond Hettinger28de64f2008-01-13 23:40:30 +0000759 def seed(self, *args, **kwds):
Raymond Hettinger23f12412004-09-13 22:23:21 +0000760 "Stub method. Not used for a system random number generator."
Raymond Hettinger356a4592004-08-30 06:14:31 +0000761 return None
Raymond Hettinger356a4592004-08-30 06:14:31 +0000762
763 def _notimplemented(self, *args, **kwds):
Raymond Hettinger23f12412004-09-13 22:23:21 +0000764 "Method should not be called for a system random number generator."
765 raise NotImplementedError('System entropy source does not have state.')
Raymond Hettinger356a4592004-08-30 06:14:31 +0000766 getstate = setstate = _notimplemented
767
Tim Peterscd804102001-01-25 20:25:57 +0000768## -------------------- test program --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000769
Raymond Hettinger62297132003-08-30 01:24:19 +0000770def _test_generator(n, func, args):
Tim Peters0c9886d2001-01-15 01:18:21 +0000771 import time
Guido van Rossumbe19ed72007-02-09 05:37:30 +0000772 print(n, 'times', func.__name__)
Raymond Hettingerb98154e2003-05-24 17:26:02 +0000773 total = 0.0
Tim Peters0c9886d2001-01-15 01:18:21 +0000774 sqsum = 0.0
775 smallest = 1e10
776 largest = -1e10
Victor Stinner8db5b542018-12-17 11:30:34 +0100777 t0 = time.perf_counter()
Tim Peters0c9886d2001-01-15 01:18:21 +0000778 for i in range(n):
Raymond Hettinger62297132003-08-30 01:24:19 +0000779 x = func(*args)
Raymond Hettingerb98154e2003-05-24 17:26:02 +0000780 total += x
Tim Peters0c9886d2001-01-15 01:18:21 +0000781 sqsum = sqsum + x*x
782 smallest = min(x, smallest)
783 largest = max(x, largest)
Victor Stinner8db5b542018-12-17 11:30:34 +0100784 t1 = time.perf_counter()
Guido van Rossumbe19ed72007-02-09 05:37:30 +0000785 print(round(t1-t0, 3), 'sec,', end=' ')
Raymond Hettingerb98154e2003-05-24 17:26:02 +0000786 avg = total/n
Tim Petersd7b5e882001-01-25 03:36:26 +0000787 stddev = _sqrt(sqsum/n - avg*avg)
Raymond Hettinger1f548142014-05-19 20:21:43 +0100788 print('avg %g, stddev %g, min %g, max %g\n' % \
Guido van Rossumbe19ed72007-02-09 05:37:30 +0000789 (avg, stddev, smallest, largest))
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000790
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000791
792def _test(N=2000):
Raymond Hettinger62297132003-08-30 01:24:19 +0000793 _test_generator(N, random, ())
794 _test_generator(N, normalvariate, (0.0, 1.0))
795 _test_generator(N, lognormvariate, (0.0, 1.0))
796 _test_generator(N, vonmisesvariate, (0.0, 1.0))
797 _test_generator(N, gammavariate, (0.01, 1.0))
798 _test_generator(N, gammavariate, (0.1, 1.0))
799 _test_generator(N, gammavariate, (0.1, 2.0))
800 _test_generator(N, gammavariate, (0.5, 1.0))
801 _test_generator(N, gammavariate, (0.9, 1.0))
802 _test_generator(N, gammavariate, (1.0, 1.0))
803 _test_generator(N, gammavariate, (2.0, 1.0))
804 _test_generator(N, gammavariate, (20.0, 1.0))
805 _test_generator(N, gammavariate, (200.0, 1.0))
806 _test_generator(N, gauss, (0.0, 1.0))
807 _test_generator(N, betavariate, (3.0, 3.0))
Christian Heimesfe337bf2008-03-23 21:54:12 +0000808 _test_generator(N, triangular, (0.0, 1.0, 1.0/3.0))
Tim Peterscd804102001-01-25 20:25:57 +0000809
Tim Peters715c4c42001-01-26 22:56:56 +0000810# Create one instance, seeded from current time, and export its methods
Raymond Hettinger40f62172002-12-29 23:03:38 +0000811# 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
Tim Petersd7b5e882001-01-25 03:36:26 +0000816_inst = Random()
817seed = _inst.seed
818random = _inst.random
819uniform = _inst.uniform
Christian Heimesfe337bf2008-03-23 21:54:12 +0000820triangular = _inst.triangular
Tim Petersd7b5e882001-01-25 03:36:26 +0000821randint = _inst.randint
822choice = _inst.choice
823randrange = _inst.randrange
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000824sample = _inst.sample
Tim Petersd7b5e882001-01-25 03:36:26 +0000825shuffle = _inst.shuffle
Raymond Hettinger28aa4a02016-09-07 00:08:44 -0700826choices = _inst.choices
Tim Petersd7b5e882001-01-25 03:36:26 +0000827normalvariate = _inst.normalvariate
828lognormvariate = _inst.lognormvariate
Tim Petersd7b5e882001-01-25 03:36:26 +0000829expovariate = _inst.expovariate
830vonmisesvariate = _inst.vonmisesvariate
831gammavariate = _inst.gammavariate
Tim Petersd7b5e882001-01-25 03:36:26 +0000832gauss = _inst.gauss
833betavariate = _inst.betavariate
834paretovariate = _inst.paretovariate
835weibullvariate = _inst.weibullvariate
836getstate = _inst.getstate
837setstate = _inst.setstate
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000838getrandbits = _inst.getrandbits
Victor Stinner9f5fe792020-04-17 19:05:35 +0200839randbytes = _inst.randbytes
Tim Petersd7b5e882001-01-25 03:36:26 +0000840
Antoine Pitrou346cbd32017-05-27 17:50:54 +0200841if hasattr(_os, "fork"):
Gregory P. Smith163468a2017-05-29 10:03:41 -0700842 _os.register_at_fork(after_in_child=_inst.seed)
Antoine Pitrou346cbd32017-05-27 17:50:54 +0200843
844
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000845if __name__ == '__main__':
Tim Petersd7b5e882001-01-25 03:36:26 +0000846 _test()