blob: 636d1024782d063da3b1c7d3f148c50652836a29 [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
Tim Petersd7b5e882001-01-25 03:36:26 +000011 generate random permutation
12
Guido van Rossume7b146f2000-02-04 15:28:42 +000013 distributions on the real line:
14 ------------------------------
Tim Petersd7b5e882001-01-25 03:36:26 +000015 uniform
Christian Heimesfe337bf2008-03-23 21:54:12 +000016 triangular
Guido van Rossume7b146f2000-02-04 15:28:42 +000017 normal (Gaussian)
18 lognormal
19 negative exponential
20 gamma
21 beta
Raymond Hettinger40f62172002-12-29 23:03:38 +000022 pareto
23 Weibull
Guido van Rossumff03b1a1994-03-09 12:55:02 +000024
Guido van Rossume7b146f2000-02-04 15:28:42 +000025 distributions on the circle (angles 0 to 2pi)
26 ---------------------------------------------
27 circular uniform
28 von Mises
29
Raymond Hettinger40f62172002-12-29 23:03:38 +000030General notes on the underlying Mersenne Twister core generator:
Guido van Rossume7b146f2000-02-04 15:28:42 +000031
Raymond Hettinger40f62172002-12-29 23:03:38 +000032* The period is 2**19937-1.
Thomas Wouters0e3f5912006-08-11 14:57:12 +000033* It is one of the most extensively tested generators in existence.
Thomas Wouters0e3f5912006-08-11 14:57:12 +000034* The random() method is implemented in C, executes in a single Python step,
35 and is, therefore, threadsafe.
Tim Peterse360d952001-01-26 10:00:39 +000036
Guido van Rossume7b146f2000-02-04 15:28:42 +000037"""
Guido van Rossumd03e1191998-05-29 17:51:31 +000038
Christian Heimesfe337bf2008-03-23 21:54:12 +000039from __future__ import division
Raymond Hettinger2f726e92003-10-05 09:09:15 +000040from warnings import warn as _warn
41from types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType
Raymond Hettinger91e27c22005-08-19 01:36:35 +000042from math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil
Tim Petersd7b5e882001-01-25 03:36:26 +000043from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
Raymond Hettingerc1c43ca2004-09-05 00:00:42 +000044from os import urandom as _urandom
45from binascii import hexlify as _hexlify
Raymond Hettinger886687d2009-02-24 11:27:15 +000046import collections as _collections
Guido van Rossumff03b1a1994-03-09 12:55:02 +000047
Raymond Hettingerf24eb352002-11-12 17:41:57 +000048__all__ = ["Random","seed","random","uniform","randint","choice","sample",
Skip Montanaro0de65802001-02-15 22:15:14 +000049 "randrange","shuffle","normalvariate","lognormvariate",
Christian Heimesfe337bf2008-03-23 21:54:12 +000050 "expovariate","vonmisesvariate","gammavariate","triangular",
Raymond Hettingerf8a52d32003-08-05 12:23:19 +000051 "gauss","betavariate","paretovariate","weibullvariate",
Raymond Hettinger28de64f2008-01-13 23:40:30 +000052 "getstate","setstate", "getrandbits",
Raymond Hettinger23f12412004-09-13 22:23:21 +000053 "SystemRandom"]
Tim Petersd7b5e882001-01-25 03:36:26 +000054
55NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
Tim Petersd7b5e882001-01-25 03:36:26 +000056TWOPI = 2.0*_pi
Tim Petersd7b5e882001-01-25 03:36:26 +000057LOG4 = _log(4.0)
Tim Petersd7b5e882001-01-25 03:36:26 +000058SG_MAGICCONST = 1.0 + _log(4.5)
Raymond Hettinger2f726e92003-10-05 09:09:15 +000059BPF = 53 # Number of bits in a float
Tim Peters7c2a85b2004-08-31 02:19:55 +000060RECIP_BPF = 2**-BPF
Tim Petersd7b5e882001-01-25 03:36:26 +000061
Raymond Hettinger356a4592004-08-30 06:14:31 +000062
Tim Petersd7b5e882001-01-25 03:36:26 +000063# Translated by Guido van Rossum from C source provided by
Raymond Hettinger40f62172002-12-29 23:03:38 +000064# Adrian Baddeley. Adapted by Raymond Hettinger for use with
Raymond Hettinger3fa19d72004-08-31 01:05:15 +000065# the Mersenne Twister and os.urandom() core generators.
Tim Petersd7b5e882001-01-25 03:36:26 +000066
Raymond Hettinger145a4a02003-01-07 10:25:55 +000067import _random
Raymond Hettinger40f62172002-12-29 23:03:38 +000068
Raymond Hettinger145a4a02003-01-07 10:25:55 +000069class Random(_random.Random):
Raymond Hettingerc32f0332002-05-23 19:44:49 +000070 """Random number generator base class used by bound module functions.
71
72 Used to instantiate instances of Random to get generators that don't
Raymond Hettinger28de64f2008-01-13 23:40:30 +000073 share state.
Raymond Hettingerc32f0332002-05-23 19:44:49 +000074
75 Class Random can also be subclassed if you want to use a different basic
76 generator of your own devising: in that case, override the following
Raymond Hettinger28de64f2008-01-13 23:40:30 +000077 methods: random(), seed(), getstate(), and setstate().
Benjamin Petersond18de0e2008-07-31 20:21:46 +000078 Optionally, implement a getrandbits() method so that randrange()
Raymond Hettinger2f726e92003-10-05 09:09:15 +000079 can cover arbitrarily large ranges.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +000080
Raymond Hettingerc32f0332002-05-23 19:44:49 +000081 """
Tim Petersd7b5e882001-01-25 03:36:26 +000082
Christian Heimescbf3b5c2007-12-03 21:02:03 +000083 VERSION = 3 # used by getstate/setstate
Tim Petersd7b5e882001-01-25 03:36:26 +000084
85 def __init__(self, x=None):
86 """Initialize an instance.
87
88 Optional argument x controls seeding, as for Random.seed().
89 """
90
91 self.seed(x)
Raymond Hettinger40f62172002-12-29 23:03:38 +000092 self.gauss_next = None
Tim Petersd7b5e882001-01-25 03:36:26 +000093
Raymond Hettingerf763a722010-09-07 00:38:15 +000094 def seed(self, a=None, version=2):
Tim Peters0de88fc2001-02-01 04:59:18 +000095 """Initialize internal state from hashable object.
Tim Petersd7b5e882001-01-25 03:36:26 +000096
Raymond Hettinger23f12412004-09-13 22:23:21 +000097 None or no argument seeds from current time or from an operating
98 system specific randomness source if available.
Tim Peters0de88fc2001-02-01 04:59:18 +000099
Raymond Hettingerf763a722010-09-07 00:38:15 +0000100 For version 2 (the default), all of the bits are used if a is a str,
101 bytes, or bytearray. For version 1, the hash() of a is used instead.
102
103 If a is an int, all bits are used.
104
Tim Petersd7b5e882001-01-25 03:36:26 +0000105 """
106
Raymond Hettinger3081d592003-08-09 18:30:57 +0000107 if a is None:
Raymond Hettingerc1c43ca2004-09-05 00:00:42 +0000108 try:
Guido van Rossume2a383d2007-01-15 16:59:06 +0000109 a = int(_hexlify(_urandom(16)), 16)
Raymond Hettingerc1c43ca2004-09-05 00:00:42 +0000110 except NotImplementedError:
Raymond Hettinger356a4592004-08-30 06:14:31 +0000111 import time
Guido van Rossume2a383d2007-01-15 16:59:06 +0000112 a = int(time.time() * 256) # use fractional seconds
Raymond Hettinger356a4592004-08-30 06:14:31 +0000113
Raymond Hettingerf763a722010-09-07 00:38:15 +0000114 if version == 2 and isinstance(a, (str, bytes, bytearray)):
115 if isinstance(a, str):
116 a = a.encode("utf8")
117 a = int(_hexlify(a), 16)
118
Guido van Rossumcd16bf62007-06-13 18:07:49 +0000119 super().seed(a)
Tim Peters46c04e12002-05-05 20:40:00 +0000120 self.gauss_next = None
121
Tim Peterscd804102001-01-25 20:25:57 +0000122 def getstate(self):
123 """Return internal state; can be passed to setstate() later."""
Guido van Rossumcd16bf62007-06-13 18:07:49 +0000124 return self.VERSION, super().getstate(), self.gauss_next
Tim Peterscd804102001-01-25 20:25:57 +0000125
126 def setstate(self, state):
127 """Restore internal state from object returned by getstate()."""
128 version = state[0]
Christian Heimescbf3b5c2007-12-03 21:02:03 +0000129 if version == 3:
Raymond Hettinger40f62172002-12-29 23:03:38 +0000130 version, internalstate, self.gauss_next = state
Guido van Rossumcd16bf62007-06-13 18:07:49 +0000131 super().setstate(internalstate)
Christian Heimescbf3b5c2007-12-03 21:02:03 +0000132 elif version == 2:
133 version, internalstate, self.gauss_next = state
134 # In version 2, the state was saved as signed ints, which causes
135 # inconsistencies between 32/64-bit systems. The state is
136 # really unsigned 32-bit ints, so we convert negative ints from
137 # version 2 to positive longs for version 3.
138 try:
139 internalstate = tuple( x % (2**32) for x in internalstate )
140 except ValueError as e:
141 raise TypeError from e
142 super(Random, self).setstate(internalstate)
Tim Peterscd804102001-01-25 20:25:57 +0000143 else:
144 raise ValueError("state with version %s passed to "
145 "Random.setstate() of version %s" %
146 (version, self.VERSION))
147
Tim Peterscd804102001-01-25 20:25:57 +0000148## ---- Methods below this point do not need to be overridden when
149## ---- subclassing for the purpose of using a different core generator.
150
151## -------------------- pickle support -------------------
152
153 def __getstate__(self): # for pickle
154 return self.getstate()
155
156 def __setstate__(self, state): # for pickle
157 self.setstate(state)
158
Raymond Hettinger5f078ff2003-06-24 20:29:04 +0000159 def __reduce__(self):
160 return self.__class__, (), self.getstate()
161
Tim Peterscd804102001-01-25 20:25:57 +0000162## -------------------- integer methods -------------------
163
Raymond Hettinger05156612010-09-07 04:44:52 +0000164 def randrange(self, start, stop=None, step=1, int=int):
Tim Petersd7b5e882001-01-25 03:36:26 +0000165 """Choose a random item from range(start, stop[, step]).
166
167 This fixes the problem with randint() which includes the
168 endpoint; in Python this is usually not what you want.
Raymond Hettinger3051cc32010-09-07 00:48:40 +0000169
170 Do not supply the 'int' and 'maxwidth' arguments.
Tim Petersd7b5e882001-01-25 03:36:26 +0000171 """
172
173 # This code is a bit messy to make it fast for the
Tim Peters9146f272002-08-16 03:41:39 +0000174 # common case while still doing adequate error checking.
Tim Petersd7b5e882001-01-25 03:36:26 +0000175 istart = int(start)
176 if istart != start:
Collin Winterce36ad82007-08-30 01:19:48 +0000177 raise ValueError("non-integer arg 1 for randrange()")
Raymond Hettinger3051cc32010-09-07 00:48:40 +0000178 if stop is None:
Tim Petersd7b5e882001-01-25 03:36:26 +0000179 if istart > 0:
Raymond Hettinger05156612010-09-07 04:44:52 +0000180 return self._randbelow(istart)
Collin Winterce36ad82007-08-30 01:19:48 +0000181 raise ValueError("empty range for randrange()")
Tim Peters9146f272002-08-16 03:41:39 +0000182
183 # stop argument supplied.
Tim Petersd7b5e882001-01-25 03:36:26 +0000184 istop = int(stop)
185 if istop != stop:
Collin Winterce36ad82007-08-30 01:19:48 +0000186 raise ValueError("non-integer stop for randrange()")
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000187 width = istop - istart
188 if step == 1 and width > 0:
Tim Peters76ca1d42003-06-19 03:46:46 +0000189 # Note that
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000190 # int(istart + self.random()*width)
Tim Peters76ca1d42003-06-19 03:46:46 +0000191 # instead would be incorrect. For example, consider istart
192 # = -2 and istop = 0. Then the guts would be in
193 # -2.0 to 0.0 exclusive on both ends (ignoring that random()
194 # might return 0.0), and because int() truncates toward 0, the
195 # final result would be -1 or 0 (instead of -2 or -1).
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000196 # istart + int(self.random()*width)
Tim Peters76ca1d42003-06-19 03:46:46 +0000197 # would also be incorrect, for a subtler reason: the RHS
198 # can return a long, and then randrange() would also return
199 # a long, but we're supposed to return an int (for backward
200 # compatibility).
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000201
Raymond Hettinger05156612010-09-07 04:44:52 +0000202 return int(istart + self._randbelow(width))
Tim Petersd7b5e882001-01-25 03:36:26 +0000203 if step == 1:
Collin Winterce36ad82007-08-30 01:19:48 +0000204 raise ValueError("empty range for randrange() (%d,%d, %d)" % (istart, istop, width))
Tim Peters9146f272002-08-16 03:41:39 +0000205
206 # Non-unit step argument supplied.
Tim Petersd7b5e882001-01-25 03:36:26 +0000207 istep = int(step)
208 if istep != step:
Collin Winterce36ad82007-08-30 01:19:48 +0000209 raise ValueError("non-integer step for randrange()")
Tim Petersd7b5e882001-01-25 03:36:26 +0000210 if istep > 0:
Raymond Hettingerffdb8bb2004-09-27 15:29:05 +0000211 n = (width + istep - 1) // istep
Tim Petersd7b5e882001-01-25 03:36:26 +0000212 elif istep < 0:
Raymond Hettingerffdb8bb2004-09-27 15:29:05 +0000213 n = (width + istep + 1) // istep
Tim Petersd7b5e882001-01-25 03:36:26 +0000214 else:
Collin Winterce36ad82007-08-30 01:19:48 +0000215 raise ValueError("zero step for randrange()")
Tim Petersd7b5e882001-01-25 03:36:26 +0000216
217 if n <= 0:
Collin Winterce36ad82007-08-30 01:19:48 +0000218 raise ValueError("empty range for randrange()")
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000219
Raymond Hettinger05156612010-09-07 04:44:52 +0000220 return istart + istep*self._randbelow(n)
Tim Petersd7b5e882001-01-25 03:36:26 +0000221
222 def randint(self, a, b):
Tim Peterscd804102001-01-25 20:25:57 +0000223 """Return random integer in range [a, b], including both end points.
Tim Petersd7b5e882001-01-25 03:36:26 +0000224 """
225
226 return self.randrange(a, b+1)
227
Raymond Hettinger05156612010-09-07 04:44:52 +0000228 def _randbelow(self, n, int=int, _maxwidth=1<<BPF, type=type,
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000229 _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType):
230 """Return a random int in the range [0,n)
231
232 Handles the case where n has more bits than returned
233 by a single call to the underlying generator.
234 """
235
236 try:
237 getrandbits = self.getrandbits
238 except AttributeError:
239 pass
240 else:
241 # Only call self.getrandbits if the original random() builtin method
242 # has not been overridden or if a new getrandbits() was supplied.
243 # This assures that the two methods correspond.
244 if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method:
Raymond Hettinger05156612010-09-07 04:44:52 +0000245 k = n.bit_length() # don't use (n-1) here because n can be 1
246 r = getrandbits(k) # 0 <= r < 2**k
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000247 while r >= n:
248 r = getrandbits(k)
249 return r
250 if n >= _maxwidth:
251 _warn("Underlying random() generator does not supply \n"
252 "enough bits to choose from a population range this large")
253 return int(self.random() * n)
254
Tim Peterscd804102001-01-25 20:25:57 +0000255## -------------------- sequence methods -------------------
256
Tim Petersd7b5e882001-01-25 03:36:26 +0000257 def choice(self, seq):
258 """Choose a random element from a non-empty sequence."""
Raymond Hettinger05156612010-09-07 04:44:52 +0000259 return seq[self._randbelow(len(seq))] # raises IndexError if seq is empty
Tim Petersd7b5e882001-01-25 03:36:26 +0000260
261 def shuffle(self, x, random=None, int=int):
262 """x, random=random.random -> shuffle list x in place; return None.
263
264 Optional arg random is a 0-argument function returning a random
265 float in [0.0, 1.0); by default, the standard random.random.
Tim Petersd7b5e882001-01-25 03:36:26 +0000266 """
267
268 if random is None:
Raymond Hettinger05156612010-09-07 04:44:52 +0000269 for i in reversed(range(1, len(x))):
270 # pick an element in x[:i+1] with which to exchange x[i]
271 j = self._randbelow(i+1)
272 x[i], x[j] = x[j], x[i]
273 else:
274 for i in reversed(range(1, len(x))):
275 # pick an element in x[:i+1] with which to exchange x[i]
276 j = int(random() * (i+1))
277 x[i], x[j] = x[j], x[i]
Tim Petersd7b5e882001-01-25 03:36:26 +0000278
Raymond Hettingerfdbe5222003-06-13 07:01:51 +0000279 def sample(self, population, k):
Raymond Hettinger1acde192008-01-14 01:00:53 +0000280 """Chooses k unique random elements from a population sequence or set.
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000281
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000282 Returns a new list containing elements from the population while
283 leaving the original population unchanged. The resulting list is
284 in selection order so that all sub-slices will also be valid random
285 samples. This allows raffle winners (the sample) to be partitioned
286 into grand prize and second place winners (the subslices).
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000287
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000288 Members of the population need not be hashable or unique. If the
289 population contains repeats, then each occurrence is a possible
290 selection in the sample.
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000291
Guido van Rossum805365e2007-05-07 22:24:25 +0000292 To choose a sample in a range of integers, use range as an argument.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000293 This is especially fast and space efficient for sampling from a
Guido van Rossum805365e2007-05-07 22:24:25 +0000294 large population: sample(range(10000000), 60)
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000295 """
296
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000297 # Sampling without replacement entails tracking either potential
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000298 # selections (the pool) in a list or previous selections in a set.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000299
Jeremy Hylton2b55d352004-02-23 17:27:57 +0000300 # When the number of selections is small compared to the
301 # population, then tracking selections is efficient, requiring
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000302 # only a small set and an occasional reselection. For
Jeremy Hylton2b55d352004-02-23 17:27:57 +0000303 # a larger number of selections, the pool tracking method is
304 # preferred since the list takes less space than the
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000305 # set and it doesn't suffer from frequent reselections.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000306
Raymond Hettinger886687d2009-02-24 11:27:15 +0000307 if isinstance(population, _collections.Set):
Raymond Hettinger1acde192008-01-14 01:00:53 +0000308 population = tuple(population)
Raymond Hettinger886687d2009-02-24 11:27:15 +0000309 if not isinstance(population, _collections.Sequence):
310 raise TypeError("Population must be a sequence or Set. For dicts, use list(d).")
Raymond Hettinger1acde192008-01-14 01:00:53 +0000311 random = self.random
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000312 n = len(population)
313 if not 0 <= k <= n:
Raymond Hettinger1acde192008-01-14 01:00:53 +0000314 raise ValueError("Sample larger than population")
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000315 result = [None] * k
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000316 setsize = 21 # size of a small set minus size of an empty list
317 if k > 5:
Tim Peters9e34c042005-08-26 15:20:46 +0000318 setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets
Raymond Hettinger1acde192008-01-14 01:00:53 +0000319 if n <= setsize:
320 # An n-length list is smaller than a k-length set
Raymond Hettinger311f4192002-11-18 09:01:24 +0000321 pool = list(population)
Guido van Rossum805365e2007-05-07 22:24:25 +0000322 for i in range(k): # invariant: non-selected at [0,n-i)
Raymond Hettinger05156612010-09-07 04:44:52 +0000323 j = self._randbelow(n-i)
Raymond Hettinger311f4192002-11-18 09:01:24 +0000324 result[i] = pool[j]
Raymond Hettinger8b9aa8d2003-01-04 05:20:33 +0000325 pool[j] = pool[n-i-1] # move non-selected item into vacancy
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000326 else:
Raymond Hettinger1acde192008-01-14 01:00:53 +0000327 selected = set()
328 selected_add = selected.add
329 for i in range(k):
Raymond Hettinger05156612010-09-07 04:44:52 +0000330 j = self._randbelow(n)
Raymond Hettinger1acde192008-01-14 01:00:53 +0000331 while j in selected:
Raymond Hettinger05156612010-09-07 04:44:52 +0000332 j = self._randbelow(n)
Raymond Hettinger1acde192008-01-14 01:00:53 +0000333 selected_add(j)
334 result[i] = population[j]
Raymond Hettinger311f4192002-11-18 09:01:24 +0000335 return result
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000336
Tim Peterscd804102001-01-25 20:25:57 +0000337## -------------------- real-valued distributions -------------------
338
339## -------------------- uniform distribution -------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000340
341 def uniform(self, a, b):
Raymond Hettingerbe40db02009-06-11 23:12:14 +0000342 "Get a random number in the range [a, b) or [a, b] depending on rounding."
Tim Petersd7b5e882001-01-25 03:36:26 +0000343 return a + (b-a) * self.random()
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000344
Christian Heimesfe337bf2008-03-23 21:54:12 +0000345## -------------------- triangular --------------------
346
347 def triangular(self, low=0.0, high=1.0, mode=None):
348 """Triangular distribution.
349
350 Continuous distribution bounded by given lower and upper limits,
351 and having a given mode value in-between.
352
353 http://en.wikipedia.org/wiki/Triangular_distribution
354
355 """
356 u = self.random()
357 c = 0.5 if mode is None else (mode - low) / (high - low)
358 if u > c:
359 u = 1.0 - u
360 c = 1.0 - c
361 low, high = high, low
362 return low + (high - low) * (u * c) ** 0.5
363
Tim Peterscd804102001-01-25 20:25:57 +0000364## -------------------- normal distribution --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000365
Tim Petersd7b5e882001-01-25 03:36:26 +0000366 def normalvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000367 """Normal distribution.
368
369 mu is the mean, and sigma is the standard deviation.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000370
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000371 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000372 # mu = mean, sigma = standard deviation
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000373
Tim Petersd7b5e882001-01-25 03:36:26 +0000374 # Uses Kinderman and Monahan method. Reference: Kinderman,
375 # A.J. and Monahan, J.F., "Computer generation of random
376 # variables using the ratio of uniform deviates", ACM Trans
377 # Math Software, 3, (1977), pp257-260.
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000378
Tim Petersd7b5e882001-01-25 03:36:26 +0000379 random = self.random
Raymond Hettinger42406e62005-04-30 09:02:51 +0000380 while 1:
Tim Peters0c9886d2001-01-15 01:18:21 +0000381 u1 = random()
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000382 u2 = 1.0 - random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000383 z = NV_MAGICCONST*(u1-0.5)/u2
384 zz = z*z/4.0
385 if zz <= -_log(u2):
386 break
387 return mu + z*sigma
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000388
Tim Peterscd804102001-01-25 20:25:57 +0000389## -------------------- lognormal distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000390
391 def lognormvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000392 """Log normal distribution.
393
394 If you take the natural logarithm of this distribution, you'll get a
395 normal distribution with mean mu and standard deviation sigma.
396 mu can have any value, and sigma must be greater than zero.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000397
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000398 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000399 return _exp(self.normalvariate(mu, sigma))
400
Tim Peterscd804102001-01-25 20:25:57 +0000401## -------------------- exponential distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000402
403 def expovariate(self, lambd):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000404 """Exponential distribution.
405
Mark Dickinson2f947362009-01-07 17:54:07 +0000406 lambd is 1.0 divided by the desired mean. It should be
407 nonzero. (The parameter would be called "lambda", but that is
408 a reserved word in Python.) Returned values range from 0 to
409 positive infinity if lambd is positive, and from negative
410 infinity to 0 if lambd is negative.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000411
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000412 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000413 # lambd: rate lambd = 1/mean
414 # ('lambda' is a Python reserved word)
415
416 random = self.random
Tim Peters0c9886d2001-01-15 01:18:21 +0000417 u = random()
418 while u <= 1e-7:
419 u = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000420 return -_log(u)/lambd
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000421
Tim Peterscd804102001-01-25 20:25:57 +0000422## -------------------- von Mises distribution --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000423
Tim Petersd7b5e882001-01-25 03:36:26 +0000424 def vonmisesvariate(self, mu, kappa):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000425 """Circular data distribution.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000426
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000427 mu is the mean angle, expressed in radians between 0 and 2*pi, and
428 kappa is the concentration parameter, which must be greater than or
429 equal to zero. If kappa is equal to zero, this distribution reduces
430 to a uniform random angle over the range 0 to 2*pi.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000431
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000432 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000433 # mu: mean angle (in radians between 0 and 2*pi)
434 # kappa: concentration parameter kappa (>= 0)
435 # if kappa = 0 generate uniform random angle
436
437 # Based upon an algorithm published in: Fisher, N.I.,
438 # "Statistical Analysis of Circular Data", Cambridge
439 # University Press, 1993.
440
441 # Thanks to Magnus Kessler for a correction to the
442 # implementation of step 4.
443
444 random = self.random
445 if kappa <= 1e-6:
446 return TWOPI * random()
447
448 a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa)
449 b = (a - _sqrt(2.0 * a))/(2.0 * kappa)
450 r = (1.0 + b * b)/(2.0 * b)
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000451
Raymond Hettinger42406e62005-04-30 09:02:51 +0000452 while 1:
Tim Peters0c9886d2001-01-15 01:18:21 +0000453 u1 = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000454
455 z = _cos(_pi * u1)
456 f = (1.0 + r * z)/(r + z)
457 c = kappa * (r - f)
458
459 u2 = random()
460
Raymond Hettinger42406e62005-04-30 09:02:51 +0000461 if u2 < c * (2.0 - c) or u2 <= c * _exp(1.0 - c):
Tim Peters0c9886d2001-01-15 01:18:21 +0000462 break
Tim Petersd7b5e882001-01-25 03:36:26 +0000463
464 u3 = random()
465 if u3 > 0.5:
466 theta = (mu % TWOPI) + _acos(f)
467 else:
468 theta = (mu % TWOPI) - _acos(f)
469
470 return theta
471
Tim Peterscd804102001-01-25 20:25:57 +0000472## -------------------- gamma distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000473
474 def gammavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000475 """Gamma distribution. Not the gamma function!
476
477 Conditions on the parameters are alpha > 0 and beta > 0.
478
479 """
Tim Peters8ac14952002-05-23 15:15:30 +0000480
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000481 # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
Tim Peters8ac14952002-05-23 15:15:30 +0000482
Guido van Rossum570764d2002-05-14 14:08:12 +0000483 # Warning: a few older sources define the gamma distribution in terms
484 # of alpha > -1.0
485 if alpha <= 0.0 or beta <= 0.0:
Collin Winterce36ad82007-08-30 01:19:48 +0000486 raise ValueError('gammavariate: alpha and beta must be > 0.0')
Tim Peters8ac14952002-05-23 15:15:30 +0000487
Tim Petersd7b5e882001-01-25 03:36:26 +0000488 random = self.random
Tim Petersd7b5e882001-01-25 03:36:26 +0000489 if alpha > 1.0:
490
491 # Uses R.C.H. Cheng, "The generation of Gamma
492 # variables with non-integral shape parameters",
493 # Applied Statistics, (1977), 26, No. 1, p71-74
494
Raymond Hettingerca6cdc22002-05-13 23:40:14 +0000495 ainv = _sqrt(2.0 * alpha - 1.0)
496 bbb = alpha - LOG4
497 ccc = alpha + ainv
Tim Peters8ac14952002-05-23 15:15:30 +0000498
Raymond Hettinger42406e62005-04-30 09:02:51 +0000499 while 1:
Tim Petersd7b5e882001-01-25 03:36:26 +0000500 u1 = random()
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000501 if not 1e-7 < u1 < .9999999:
502 continue
503 u2 = 1.0 - random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000504 v = _log(u1/(1.0-u1))/ainv
505 x = alpha*_exp(v)
506 z = u1*u1*u2
507 r = bbb+ccc*v-x
508 if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z):
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000509 return x * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000510
511 elif alpha == 1.0:
512 # expovariate(1)
513 u = random()
514 while u <= 1e-7:
515 u = random()
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000516 return -_log(u) * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000517
518 else: # alpha is between 0 and 1 (exclusive)
519
520 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
521
Raymond Hettinger42406e62005-04-30 09:02:51 +0000522 while 1:
Tim Petersd7b5e882001-01-25 03:36:26 +0000523 u = random()
524 b = (_e + alpha)/_e
525 p = b*u
526 if p <= 1.0:
Raymond Hettinger42406e62005-04-30 09:02:51 +0000527 x = p ** (1.0/alpha)
Tim Petersd7b5e882001-01-25 03:36:26 +0000528 else:
Tim Petersd7b5e882001-01-25 03:36:26 +0000529 x = -_log((b-p)/alpha)
530 u1 = random()
Raymond Hettinger42406e62005-04-30 09:02:51 +0000531 if p > 1.0:
532 if u1 <= x ** (alpha - 1.0):
533 break
534 elif u1 <= _exp(-x):
Tim Petersd7b5e882001-01-25 03:36:26 +0000535 break
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000536 return x * beta
537
Tim Peterscd804102001-01-25 20:25:57 +0000538## -------------------- Gauss (faster alternative) --------------------
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000539
Tim Petersd7b5e882001-01-25 03:36:26 +0000540 def gauss(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000541 """Gaussian distribution.
542
543 mu is the mean, and sigma is the standard deviation. This is
544 slightly faster than the normalvariate() function.
545
546 Not thread-safe without a lock around calls.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000547
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000548 """
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000549
Tim Petersd7b5e882001-01-25 03:36:26 +0000550 # When x and y are two variables from [0, 1), uniformly
551 # distributed, then
552 #
553 # cos(2*pi*x)*sqrt(-2*log(1-y))
554 # sin(2*pi*x)*sqrt(-2*log(1-y))
555 #
556 # are two *independent* variables with normal distribution
557 # (mu = 0, sigma = 1).
558 # (Lambert Meertens)
559 # (corrected version; bug discovered by Mike Miller, fixed by LM)
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000560
Tim Petersd7b5e882001-01-25 03:36:26 +0000561 # Multithreading note: When two threads call this function
562 # simultaneously, it is possible that they will receive the
563 # same return value. The window is very small though. To
564 # avoid this, you have to use a lock around all calls. (I
565 # didn't want to slow this down in the serial case by using a
566 # lock here.)
Guido van Rossumd03e1191998-05-29 17:51:31 +0000567
Tim Petersd7b5e882001-01-25 03:36:26 +0000568 random = self.random
569 z = self.gauss_next
570 self.gauss_next = None
571 if z is None:
572 x2pi = random() * TWOPI
573 g2rad = _sqrt(-2.0 * _log(1.0 - random()))
574 z = _cos(x2pi) * g2rad
575 self.gauss_next = _sin(x2pi) * g2rad
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000576
Tim Petersd7b5e882001-01-25 03:36:26 +0000577 return mu + z*sigma
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000578
Tim Peterscd804102001-01-25 20:25:57 +0000579## -------------------- beta --------------------
Tim Peters85e2e472001-01-26 06:49:56 +0000580## See
581## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470
582## for Ivan Frohne's insightful analysis of why the original implementation:
583##
584## def betavariate(self, alpha, beta):
585## # Discrete Event Simulation in C, pp 87-88.
586##
587## y = self.expovariate(alpha)
588## z = self.expovariate(1.0/beta)
589## return z/(y+z)
590##
591## was dead wrong, and how it probably got that way.
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000592
Tim Petersd7b5e882001-01-25 03:36:26 +0000593 def betavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000594 """Beta distribution.
595
Thomas Woutersb2137042007-02-01 18:02:27 +0000596 Conditions on the parameters are alpha > 0 and beta > 0.
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000597 Returned values range between 0 and 1.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000598
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000599 """
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000600
Tim Peters85e2e472001-01-26 06:49:56 +0000601 # This version due to Janne Sinkkonen, and matches all the std
602 # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
603 y = self.gammavariate(alpha, 1.)
604 if y == 0:
605 return 0.0
606 else:
607 return y / (y + self.gammavariate(beta, 1.))
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000608
Tim Peterscd804102001-01-25 20:25:57 +0000609## -------------------- Pareto --------------------
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000610
Tim Petersd7b5e882001-01-25 03:36:26 +0000611 def paretovariate(self, alpha):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000612 """Pareto distribution. alpha is the shape parameter."""
Tim Petersd7b5e882001-01-25 03:36:26 +0000613 # Jain, pg. 495
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000614
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000615 u = 1.0 - self.random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000616 return 1.0 / pow(u, 1.0/alpha)
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000617
Tim Peterscd804102001-01-25 20:25:57 +0000618## -------------------- Weibull --------------------
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000619
Tim Petersd7b5e882001-01-25 03:36:26 +0000620 def weibullvariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000621 """Weibull distribution.
622
623 alpha is the scale parameter and beta is the shape parameter.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000624
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000625 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000626 # Jain, pg. 499; bug fix courtesy Bill Arms
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000627
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000628 u = 1.0 - self.random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000629 return alpha * pow(-_log(u), 1.0/beta)
Guido van Rossum6c395ba1999-08-18 13:53:28 +0000630
Raymond Hettinger23f12412004-09-13 22:23:21 +0000631## --------------- Operating System Random Source ------------------
Raymond Hettinger356a4592004-08-30 06:14:31 +0000632
Raymond Hettinger23f12412004-09-13 22:23:21 +0000633class SystemRandom(Random):
634 """Alternate random number generator using sources provided
635 by the operating system (such as /dev/urandom on Unix or
636 CryptGenRandom on Windows).
Raymond Hettinger356a4592004-08-30 06:14:31 +0000637
638 Not available on all systems (see os.urandom() for details).
639 """
640
641 def random(self):
642 """Get the next random number in the range [0.0, 1.0)."""
Guido van Rossume2a383d2007-01-15 16:59:06 +0000643 return (int(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF
Raymond Hettinger356a4592004-08-30 06:14:31 +0000644
645 def getrandbits(self, k):
646 """getrandbits(k) -> x. Generates a long int with k random bits."""
Raymond Hettinger356a4592004-08-30 06:14:31 +0000647 if k <= 0:
648 raise ValueError('number of bits must be greater than zero')
649 if k != int(k):
650 raise TypeError('number of bits should be an integer')
651 bytes = (k + 7) // 8 # bits / 8 and rounded up
Guido van Rossume2a383d2007-01-15 16:59:06 +0000652 x = int(_hexlify(_urandom(bytes)), 16)
Raymond Hettinger356a4592004-08-30 06:14:31 +0000653 return x >> (bytes * 8 - k) # trim excess bits
654
Raymond Hettinger28de64f2008-01-13 23:40:30 +0000655 def seed(self, *args, **kwds):
Raymond Hettinger23f12412004-09-13 22:23:21 +0000656 "Stub method. Not used for a system random number generator."
Raymond Hettinger356a4592004-08-30 06:14:31 +0000657 return None
Raymond Hettinger356a4592004-08-30 06:14:31 +0000658
659 def _notimplemented(self, *args, **kwds):
Raymond Hettinger23f12412004-09-13 22:23:21 +0000660 "Method should not be called for a system random number generator."
661 raise NotImplementedError('System entropy source does not have state.')
Raymond Hettinger356a4592004-08-30 06:14:31 +0000662 getstate = setstate = _notimplemented
663
Tim Peterscd804102001-01-25 20:25:57 +0000664## -------------------- test program --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000665
Raymond Hettinger62297132003-08-30 01:24:19 +0000666def _test_generator(n, func, args):
Tim Peters0c9886d2001-01-15 01:18:21 +0000667 import time
Guido van Rossumbe19ed72007-02-09 05:37:30 +0000668 print(n, 'times', func.__name__)
Raymond Hettingerb98154e2003-05-24 17:26:02 +0000669 total = 0.0
Tim Peters0c9886d2001-01-15 01:18:21 +0000670 sqsum = 0.0
671 smallest = 1e10
672 largest = -1e10
673 t0 = time.time()
674 for i in range(n):
Raymond Hettinger62297132003-08-30 01:24:19 +0000675 x = func(*args)
Raymond Hettingerb98154e2003-05-24 17:26:02 +0000676 total += x
Tim Peters0c9886d2001-01-15 01:18:21 +0000677 sqsum = sqsum + x*x
678 smallest = min(x, smallest)
679 largest = max(x, largest)
680 t1 = time.time()
Guido van Rossumbe19ed72007-02-09 05:37:30 +0000681 print(round(t1-t0, 3), 'sec,', end=' ')
Raymond Hettingerb98154e2003-05-24 17:26:02 +0000682 avg = total/n
Tim Petersd7b5e882001-01-25 03:36:26 +0000683 stddev = _sqrt(sqsum/n - avg*avg)
Guido van Rossumbe19ed72007-02-09 05:37:30 +0000684 print('avg %g, stddev %g, min %g, max %g' % \
685 (avg, stddev, smallest, largest))
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000686
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000687
688def _test(N=2000):
Raymond Hettinger62297132003-08-30 01:24:19 +0000689 _test_generator(N, random, ())
690 _test_generator(N, normalvariate, (0.0, 1.0))
691 _test_generator(N, lognormvariate, (0.0, 1.0))
692 _test_generator(N, vonmisesvariate, (0.0, 1.0))
693 _test_generator(N, gammavariate, (0.01, 1.0))
694 _test_generator(N, gammavariate, (0.1, 1.0))
695 _test_generator(N, gammavariate, (0.1, 2.0))
696 _test_generator(N, gammavariate, (0.5, 1.0))
697 _test_generator(N, gammavariate, (0.9, 1.0))
698 _test_generator(N, gammavariate, (1.0, 1.0))
699 _test_generator(N, gammavariate, (2.0, 1.0))
700 _test_generator(N, gammavariate, (20.0, 1.0))
701 _test_generator(N, gammavariate, (200.0, 1.0))
702 _test_generator(N, gauss, (0.0, 1.0))
703 _test_generator(N, betavariate, (3.0, 3.0))
Christian Heimesfe337bf2008-03-23 21:54:12 +0000704 _test_generator(N, triangular, (0.0, 1.0, 1.0/3.0))
Tim Peterscd804102001-01-25 20:25:57 +0000705
Tim Peters715c4c42001-01-26 22:56:56 +0000706# Create one instance, seeded from current time, and export its methods
Raymond Hettinger40f62172002-12-29 23:03:38 +0000707# as module-level functions. The functions share state across all uses
708#(both in the user's code and in the Python libraries), but that's fine
709# for most programs and is easier for the casual user than making them
710# instantiate their own Random() instance.
711
Tim Petersd7b5e882001-01-25 03:36:26 +0000712_inst = Random()
713seed = _inst.seed
714random = _inst.random
715uniform = _inst.uniform
Christian Heimesfe337bf2008-03-23 21:54:12 +0000716triangular = _inst.triangular
Tim Petersd7b5e882001-01-25 03:36:26 +0000717randint = _inst.randint
718choice = _inst.choice
719randrange = _inst.randrange
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000720sample = _inst.sample
Tim Petersd7b5e882001-01-25 03:36:26 +0000721shuffle = _inst.shuffle
722normalvariate = _inst.normalvariate
723lognormvariate = _inst.lognormvariate
Tim Petersd7b5e882001-01-25 03:36:26 +0000724expovariate = _inst.expovariate
725vonmisesvariate = _inst.vonmisesvariate
726gammavariate = _inst.gammavariate
Tim Petersd7b5e882001-01-25 03:36:26 +0000727gauss = _inst.gauss
728betavariate = _inst.betavariate
729paretovariate = _inst.paretovariate
730weibullvariate = _inst.weibullvariate
731getstate = _inst.getstate
732setstate = _inst.setstate
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000733getrandbits = _inst.getrandbits
Tim Petersd7b5e882001-01-25 03:36:26 +0000734
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000735if __name__ == '__main__':
Tim Petersd7b5e882001-01-25 03:36:26 +0000736 _test()