blob: 592e4b899a17b9d218ec0d9f058785951389aeb7 [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
Tim Peters0de88fc2001-02-01 04:59:18 +000094 def seed(self, a=None):
95 """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
Mark Dickinson5c2db372009-12-05 20:28:34 +0000100 If a is not None or an int, hash(a) is used instead.
Tim Petersd7b5e882001-01-25 03:36:26 +0000101 """
102
Raymond Hettinger3081d592003-08-09 18:30:57 +0000103 if a is None:
Raymond Hettingerc1c43ca2004-09-05 00:00:42 +0000104 try:
Guido van Rossume2a383d2007-01-15 16:59:06 +0000105 a = int(_hexlify(_urandom(16)), 16)
Raymond Hettingerc1c43ca2004-09-05 00:00:42 +0000106 except NotImplementedError:
Raymond Hettinger356a4592004-08-30 06:14:31 +0000107 import time
Guido van Rossume2a383d2007-01-15 16:59:06 +0000108 a = int(time.time() * 256) # use fractional seconds
Raymond Hettinger356a4592004-08-30 06:14:31 +0000109
Guido van Rossumcd16bf62007-06-13 18:07:49 +0000110 super().seed(a)
Tim Peters46c04e12002-05-05 20:40:00 +0000111 self.gauss_next = None
112
Tim Peterscd804102001-01-25 20:25:57 +0000113 def getstate(self):
114 """Return internal state; can be passed to setstate() later."""
Guido van Rossumcd16bf62007-06-13 18:07:49 +0000115 return self.VERSION, super().getstate(), self.gauss_next
Tim Peterscd804102001-01-25 20:25:57 +0000116
117 def setstate(self, state):
118 """Restore internal state from object returned by getstate()."""
119 version = state[0]
Christian Heimescbf3b5c2007-12-03 21:02:03 +0000120 if version == 3:
Raymond Hettinger40f62172002-12-29 23:03:38 +0000121 version, internalstate, self.gauss_next = state
Guido van Rossumcd16bf62007-06-13 18:07:49 +0000122 super().setstate(internalstate)
Christian Heimescbf3b5c2007-12-03 21:02:03 +0000123 elif version == 2:
124 version, internalstate, self.gauss_next = state
125 # In version 2, the state was saved as signed ints, which causes
126 # inconsistencies between 32/64-bit systems. The state is
127 # really unsigned 32-bit ints, so we convert negative ints from
128 # version 2 to positive longs for version 3.
129 try:
130 internalstate = tuple( x % (2**32) for x in internalstate )
131 except ValueError as e:
132 raise TypeError from e
133 super(Random, self).setstate(internalstate)
Tim Peterscd804102001-01-25 20:25:57 +0000134 else:
135 raise ValueError("state with version %s passed to "
136 "Random.setstate() of version %s" %
137 (version, self.VERSION))
138
Tim Peterscd804102001-01-25 20:25:57 +0000139## ---- Methods below this point do not need to be overridden when
140## ---- subclassing for the purpose of using a different core generator.
141
142## -------------------- pickle support -------------------
143
144 def __getstate__(self): # for pickle
145 return self.getstate()
146
147 def __setstate__(self, state): # for pickle
148 self.setstate(state)
149
Raymond Hettinger5f078ff2003-06-24 20:29:04 +0000150 def __reduce__(self):
151 return self.__class__, (), self.getstate()
152
Tim Peterscd804102001-01-25 20:25:57 +0000153## -------------------- integer methods -------------------
154
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000155 def randrange(self, start, stop=None, step=1, int=int, default=None,
Guido van Rossume2a383d2007-01-15 16:59:06 +0000156 maxwidth=1<<BPF):
Tim Petersd7b5e882001-01-25 03:36:26 +0000157 """Choose a random item from range(start, stop[, step]).
158
159 This fixes the problem with randint() which includes the
160 endpoint; in Python this is usually not what you want.
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000161 Do not supply the 'int', 'default', and 'maxwidth' arguments.
Tim Petersd7b5e882001-01-25 03:36:26 +0000162 """
163
164 # This code is a bit messy to make it fast for the
Tim Peters9146f272002-08-16 03:41:39 +0000165 # common case while still doing adequate error checking.
Tim Petersd7b5e882001-01-25 03:36:26 +0000166 istart = int(start)
167 if istart != start:
Collin Winterce36ad82007-08-30 01:19:48 +0000168 raise ValueError("non-integer arg 1 for randrange()")
Tim Petersd7b5e882001-01-25 03:36:26 +0000169 if stop is default:
170 if istart > 0:
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000171 if istart >= maxwidth:
172 return self._randbelow(istart)
Tim Petersd7b5e882001-01-25 03:36:26 +0000173 return int(self.random() * istart)
Collin Winterce36ad82007-08-30 01:19:48 +0000174 raise ValueError("empty range for randrange()")
Tim Peters9146f272002-08-16 03:41:39 +0000175
176 # stop argument supplied.
Tim Petersd7b5e882001-01-25 03:36:26 +0000177 istop = int(stop)
178 if istop != stop:
Collin Winterce36ad82007-08-30 01:19:48 +0000179 raise ValueError("non-integer stop for randrange()")
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000180 width = istop - istart
181 if step == 1 and width > 0:
Tim Peters76ca1d42003-06-19 03:46:46 +0000182 # Note that
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000183 # int(istart + self.random()*width)
Tim Peters76ca1d42003-06-19 03:46:46 +0000184 # instead would be incorrect. For example, consider istart
185 # = -2 and istop = 0. Then the guts would be in
186 # -2.0 to 0.0 exclusive on both ends (ignoring that random()
187 # might return 0.0), and because int() truncates toward 0, the
188 # final result would be -1 or 0 (instead of -2 or -1).
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000189 # istart + int(self.random()*width)
Tim Peters76ca1d42003-06-19 03:46:46 +0000190 # would also be incorrect, for a subtler reason: the RHS
191 # can return a long, and then randrange() would also return
192 # a long, but we're supposed to return an int (for backward
193 # compatibility).
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000194
195 if width >= maxwidth:
Tim Peters58eb11c2004-01-18 20:29:55 +0000196 return int(istart + self._randbelow(width))
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000197 return int(istart + int(self.random()*width))
Tim Petersd7b5e882001-01-25 03:36:26 +0000198 if step == 1:
Collin Winterce36ad82007-08-30 01:19:48 +0000199 raise ValueError("empty range for randrange() (%d,%d, %d)" % (istart, istop, width))
Tim Peters9146f272002-08-16 03:41:39 +0000200
201 # Non-unit step argument supplied.
Tim Petersd7b5e882001-01-25 03:36:26 +0000202 istep = int(step)
203 if istep != step:
Collin Winterce36ad82007-08-30 01:19:48 +0000204 raise ValueError("non-integer step for randrange()")
Tim Petersd7b5e882001-01-25 03:36:26 +0000205 if istep > 0:
Raymond Hettingerffdb8bb2004-09-27 15:29:05 +0000206 n = (width + istep - 1) // istep
Tim Petersd7b5e882001-01-25 03:36:26 +0000207 elif istep < 0:
Raymond Hettingerffdb8bb2004-09-27 15:29:05 +0000208 n = (width + istep + 1) // istep
Tim Petersd7b5e882001-01-25 03:36:26 +0000209 else:
Collin Winterce36ad82007-08-30 01:19:48 +0000210 raise ValueError("zero step for randrange()")
Tim Petersd7b5e882001-01-25 03:36:26 +0000211
212 if n <= 0:
Collin Winterce36ad82007-08-30 01:19:48 +0000213 raise ValueError("empty range for randrange()")
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000214
215 if n >= maxwidth:
Thomas Wouters902d6eb2007-01-09 23:18:33 +0000216 return istart + istep*self._randbelow(n)
Tim Petersd7b5e882001-01-25 03:36:26 +0000217 return istart + istep*int(self.random() * n)
218
219 def randint(self, a, b):
Tim Peterscd804102001-01-25 20:25:57 +0000220 """Return random integer in range [a, b], including both end points.
Tim Petersd7b5e882001-01-25 03:36:26 +0000221 """
222
223 return self.randrange(a, b+1)
224
Guido van Rossume2a383d2007-01-15 16:59:06 +0000225 def _randbelow(self, n, _log=_log, int=int, _maxwidth=1<<BPF,
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000226 _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType):
227 """Return a random int in the range [0,n)
228
229 Handles the case where n has more bits than returned
230 by a single call to the underlying generator.
231 """
232
233 try:
234 getrandbits = self.getrandbits
235 except AttributeError:
236 pass
237 else:
238 # Only call self.getrandbits if the original random() builtin method
239 # has not been overridden or if a new getrandbits() was supplied.
240 # This assures that the two methods correspond.
241 if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method:
242 k = int(1.00001 + _log(n-1, 2.0)) # 2**k > n-1 > 2**(k-2)
243 r = getrandbits(k)
244 while r >= n:
245 r = getrandbits(k)
246 return r
247 if n >= _maxwidth:
248 _warn("Underlying random() generator does not supply \n"
249 "enough bits to choose from a population range this large")
250 return int(self.random() * n)
251
Tim Peterscd804102001-01-25 20:25:57 +0000252## -------------------- sequence methods -------------------
253
Tim Petersd7b5e882001-01-25 03:36:26 +0000254 def choice(self, seq):
255 """Choose a random element from a non-empty sequence."""
Raymond Hettinger5dae5052004-06-07 02:07:15 +0000256 return seq[int(self.random() * len(seq))] # raises IndexError if seq is empty
Tim Petersd7b5e882001-01-25 03:36:26 +0000257
258 def shuffle(self, x, random=None, int=int):
259 """x, random=random.random -> shuffle list x in place; return None.
260
261 Optional arg random is a 0-argument function returning a random
262 float in [0.0, 1.0); by default, the standard random.random.
Tim Petersd7b5e882001-01-25 03:36:26 +0000263 """
264
265 if random is None:
266 random = self.random
Guido van Rossum805365e2007-05-07 22:24:25 +0000267 for i in reversed(range(1, len(x))):
Tim Peterscd804102001-01-25 20:25:57 +0000268 # pick an element in x[:i+1] with which to exchange x[i]
Tim Petersd7b5e882001-01-25 03:36:26 +0000269 j = int(random() * (i+1))
270 x[i], x[j] = x[j], x[i]
271
Raymond Hettingerfdbe5222003-06-13 07:01:51 +0000272 def sample(self, population, k):
Raymond Hettinger1acde192008-01-14 01:00:53 +0000273 """Chooses k unique random elements from a population sequence or set.
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000274
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000275 Returns a new list containing elements from the population while
276 leaving the original population unchanged. The resulting list is
277 in selection order so that all sub-slices will also be valid random
278 samples. This allows raffle winners (the sample) to be partitioned
279 into grand prize and second place winners (the subslices).
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000280
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000281 Members of the population need not be hashable or unique. If the
282 population contains repeats, then each occurrence is a possible
283 selection in the sample.
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000284
Guido van Rossum805365e2007-05-07 22:24:25 +0000285 To choose a sample in a range of integers, use range as an argument.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000286 This is especially fast and space efficient for sampling from a
Guido van Rossum805365e2007-05-07 22:24:25 +0000287 large population: sample(range(10000000), 60)
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000288 """
289
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000290 # Sampling without replacement entails tracking either potential
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000291 # selections (the pool) in a list or previous selections in a set.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000292
Jeremy Hylton2b55d352004-02-23 17:27:57 +0000293 # When the number of selections is small compared to the
294 # population, then tracking selections is efficient, requiring
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000295 # only a small set and an occasional reselection. For
Jeremy Hylton2b55d352004-02-23 17:27:57 +0000296 # a larger number of selections, the pool tracking method is
297 # preferred since the list takes less space than the
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000298 # set and it doesn't suffer from frequent reselections.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000299
Raymond Hettinger886687d2009-02-24 11:27:15 +0000300 if isinstance(population, _collections.Set):
Raymond Hettinger1acde192008-01-14 01:00:53 +0000301 population = tuple(population)
Raymond Hettinger886687d2009-02-24 11:27:15 +0000302 if not isinstance(population, _collections.Sequence):
303 raise TypeError("Population must be a sequence or Set. For dicts, use list(d).")
Raymond Hettinger1acde192008-01-14 01:00:53 +0000304 random = self.random
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000305 n = len(population)
306 if not 0 <= k <= n:
Raymond Hettinger1acde192008-01-14 01:00:53 +0000307 raise ValueError("Sample larger than population")
Raymond Hettingerfdbe5222003-06-13 07:01:51 +0000308 _int = int
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000309 result = [None] * k
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000310 setsize = 21 # size of a small set minus size of an empty list
311 if k > 5:
Tim Peters9e34c042005-08-26 15:20:46 +0000312 setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets
Raymond Hettinger1acde192008-01-14 01:00:53 +0000313 if n <= setsize:
314 # An n-length list is smaller than a k-length set
Raymond Hettinger311f4192002-11-18 09:01:24 +0000315 pool = list(population)
Guido van Rossum805365e2007-05-07 22:24:25 +0000316 for i in range(k): # invariant: non-selected at [0,n-i)
Raymond Hettingerfdbe5222003-06-13 07:01:51 +0000317 j = _int(random() * (n-i))
Raymond Hettinger311f4192002-11-18 09:01:24 +0000318 result[i] = pool[j]
Raymond Hettinger8b9aa8d2003-01-04 05:20:33 +0000319 pool[j] = pool[n-i-1] # move non-selected item into vacancy
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000320 else:
Raymond Hettinger1acde192008-01-14 01:00:53 +0000321 selected = set()
322 selected_add = selected.add
323 for i in range(k):
324 j = _int(random() * n)
325 while j in selected:
Raymond Hettingerfdbe5222003-06-13 07:01:51 +0000326 j = _int(random() * n)
Raymond Hettinger1acde192008-01-14 01:00:53 +0000327 selected_add(j)
328 result[i] = population[j]
Raymond Hettinger311f4192002-11-18 09:01:24 +0000329 return result
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000330
Tim Peterscd804102001-01-25 20:25:57 +0000331## -------------------- real-valued distributions -------------------
332
333## -------------------- uniform distribution -------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000334
335 def uniform(self, a, b):
Raymond Hettingerbe40db02009-06-11 23:12:14 +0000336 "Get a random number in the range [a, b) or [a, b] depending on rounding."
Tim Petersd7b5e882001-01-25 03:36:26 +0000337 return a + (b-a) * self.random()
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000338
Christian Heimesfe337bf2008-03-23 21:54:12 +0000339## -------------------- triangular --------------------
340
341 def triangular(self, low=0.0, high=1.0, mode=None):
342 """Triangular distribution.
343
344 Continuous distribution bounded by given lower and upper limits,
345 and having a given mode value in-between.
346
347 http://en.wikipedia.org/wiki/Triangular_distribution
348
349 """
350 u = self.random()
351 c = 0.5 if mode is None else (mode - low) / (high - low)
352 if u > c:
353 u = 1.0 - u
354 c = 1.0 - c
355 low, high = high, low
356 return low + (high - low) * (u * c) ** 0.5
357
Tim Peterscd804102001-01-25 20:25:57 +0000358## -------------------- normal distribution --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000359
Tim Petersd7b5e882001-01-25 03:36:26 +0000360 def normalvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000361 """Normal distribution.
362
363 mu is the mean, and sigma is the standard deviation.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000364
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000365 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000366 # mu = mean, sigma = standard deviation
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000367
Tim Petersd7b5e882001-01-25 03:36:26 +0000368 # Uses Kinderman and Monahan method. Reference: Kinderman,
369 # A.J. and Monahan, J.F., "Computer generation of random
370 # variables using the ratio of uniform deviates", ACM Trans
371 # Math Software, 3, (1977), pp257-260.
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000372
Tim Petersd7b5e882001-01-25 03:36:26 +0000373 random = self.random
Raymond Hettinger42406e62005-04-30 09:02:51 +0000374 while 1:
Tim Peters0c9886d2001-01-15 01:18:21 +0000375 u1 = random()
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000376 u2 = 1.0 - random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000377 z = NV_MAGICCONST*(u1-0.5)/u2
378 zz = z*z/4.0
379 if zz <= -_log(u2):
380 break
381 return mu + z*sigma
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000382
Tim Peterscd804102001-01-25 20:25:57 +0000383## -------------------- lognormal distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000384
385 def lognormvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000386 """Log normal distribution.
387
388 If you take the natural logarithm of this distribution, you'll get a
389 normal distribution with mean mu and standard deviation sigma.
390 mu can have any value, and sigma must be greater than zero.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000391
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000392 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000393 return _exp(self.normalvariate(mu, sigma))
394
Tim Peterscd804102001-01-25 20:25:57 +0000395## -------------------- exponential distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000396
397 def expovariate(self, lambd):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000398 """Exponential distribution.
399
Mark Dickinson2f947362009-01-07 17:54:07 +0000400 lambd is 1.0 divided by the desired mean. It should be
401 nonzero. (The parameter would be called "lambda", but that is
402 a reserved word in Python.) Returned values range from 0 to
403 positive infinity if lambd is positive, and from negative
404 infinity to 0 if lambd is negative.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000405
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000406 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000407 # lambd: rate lambd = 1/mean
408 # ('lambda' is a Python reserved word)
409
410 random = self.random
Tim Peters0c9886d2001-01-15 01:18:21 +0000411 u = random()
412 while u <= 1e-7:
413 u = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000414 return -_log(u)/lambd
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000415
Tim Peterscd804102001-01-25 20:25:57 +0000416## -------------------- von Mises distribution --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000417
Tim Petersd7b5e882001-01-25 03:36:26 +0000418 def vonmisesvariate(self, mu, kappa):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000419 """Circular data distribution.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000420
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000421 mu is the mean angle, expressed in radians between 0 and 2*pi, and
422 kappa is the concentration parameter, which must be greater than or
423 equal to zero. If kappa is equal to zero, this distribution reduces
424 to a uniform random angle over the range 0 to 2*pi.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000425
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000426 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000427 # mu: mean angle (in radians between 0 and 2*pi)
428 # kappa: concentration parameter kappa (>= 0)
429 # if kappa = 0 generate uniform random angle
430
431 # Based upon an algorithm published in: Fisher, N.I.,
432 # "Statistical Analysis of Circular Data", Cambridge
433 # University Press, 1993.
434
435 # Thanks to Magnus Kessler for a correction to the
436 # implementation of step 4.
437
438 random = self.random
439 if kappa <= 1e-6:
440 return TWOPI * random()
441
442 a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa)
443 b = (a - _sqrt(2.0 * a))/(2.0 * kappa)
444 r = (1.0 + b * b)/(2.0 * b)
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000445
Raymond Hettinger42406e62005-04-30 09:02:51 +0000446 while 1:
Tim Peters0c9886d2001-01-15 01:18:21 +0000447 u1 = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000448
449 z = _cos(_pi * u1)
450 f = (1.0 + r * z)/(r + z)
451 c = kappa * (r - f)
452
453 u2 = random()
454
Raymond Hettinger42406e62005-04-30 09:02:51 +0000455 if u2 < c * (2.0 - c) or u2 <= c * _exp(1.0 - c):
Tim Peters0c9886d2001-01-15 01:18:21 +0000456 break
Tim Petersd7b5e882001-01-25 03:36:26 +0000457
458 u3 = random()
459 if u3 > 0.5:
460 theta = (mu % TWOPI) + _acos(f)
461 else:
462 theta = (mu % TWOPI) - _acos(f)
463
464 return theta
465
Tim Peterscd804102001-01-25 20:25:57 +0000466## -------------------- gamma distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000467
468 def gammavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000469 """Gamma distribution. Not the gamma function!
470
471 Conditions on the parameters are alpha > 0 and beta > 0.
472
473 """
Tim Peters8ac14952002-05-23 15:15:30 +0000474
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000475 # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
Tim Peters8ac14952002-05-23 15:15:30 +0000476
Guido van Rossum570764d2002-05-14 14:08:12 +0000477 # Warning: a few older sources define the gamma distribution in terms
478 # of alpha > -1.0
479 if alpha <= 0.0 or beta <= 0.0:
Collin Winterce36ad82007-08-30 01:19:48 +0000480 raise ValueError('gammavariate: alpha and beta must be > 0.0')
Tim Peters8ac14952002-05-23 15:15:30 +0000481
Tim Petersd7b5e882001-01-25 03:36:26 +0000482 random = self.random
Tim Petersd7b5e882001-01-25 03:36:26 +0000483 if alpha > 1.0:
484
485 # Uses R.C.H. Cheng, "The generation of Gamma
486 # variables with non-integral shape parameters",
487 # Applied Statistics, (1977), 26, No. 1, p71-74
488
Raymond Hettingerca6cdc22002-05-13 23:40:14 +0000489 ainv = _sqrt(2.0 * alpha - 1.0)
490 bbb = alpha - LOG4
491 ccc = alpha + ainv
Tim Peters8ac14952002-05-23 15:15:30 +0000492
Raymond Hettinger42406e62005-04-30 09:02:51 +0000493 while 1:
Tim Petersd7b5e882001-01-25 03:36:26 +0000494 u1 = random()
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000495 if not 1e-7 < u1 < .9999999:
496 continue
497 u2 = 1.0 - random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000498 v = _log(u1/(1.0-u1))/ainv
499 x = alpha*_exp(v)
500 z = u1*u1*u2
501 r = bbb+ccc*v-x
502 if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z):
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000503 return x * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000504
505 elif alpha == 1.0:
506 # expovariate(1)
507 u = random()
508 while u <= 1e-7:
509 u = random()
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000510 return -_log(u) * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000511
512 else: # alpha is between 0 and 1 (exclusive)
513
514 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
515
Raymond Hettinger42406e62005-04-30 09:02:51 +0000516 while 1:
Tim Petersd7b5e882001-01-25 03:36:26 +0000517 u = random()
518 b = (_e + alpha)/_e
519 p = b*u
520 if p <= 1.0:
Raymond Hettinger42406e62005-04-30 09:02:51 +0000521 x = p ** (1.0/alpha)
Tim Petersd7b5e882001-01-25 03:36:26 +0000522 else:
Tim Petersd7b5e882001-01-25 03:36:26 +0000523 x = -_log((b-p)/alpha)
524 u1 = random()
Raymond Hettinger42406e62005-04-30 09:02:51 +0000525 if p > 1.0:
526 if u1 <= x ** (alpha - 1.0):
527 break
528 elif u1 <= _exp(-x):
Tim Petersd7b5e882001-01-25 03:36:26 +0000529 break
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000530 return x * beta
531
Tim Peterscd804102001-01-25 20:25:57 +0000532## -------------------- Gauss (faster alternative) --------------------
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000533
Tim Petersd7b5e882001-01-25 03:36:26 +0000534 def gauss(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000535 """Gaussian distribution.
536
537 mu is the mean, and sigma is the standard deviation. This is
538 slightly faster than the normalvariate() function.
539
540 Not thread-safe without a lock around calls.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000541
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000542 """
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000543
Tim Petersd7b5e882001-01-25 03:36:26 +0000544 # When x and y are two variables from [0, 1), uniformly
545 # distributed, then
546 #
547 # cos(2*pi*x)*sqrt(-2*log(1-y))
548 # sin(2*pi*x)*sqrt(-2*log(1-y))
549 #
550 # are two *independent* variables with normal distribution
551 # (mu = 0, sigma = 1).
552 # (Lambert Meertens)
553 # (corrected version; bug discovered by Mike Miller, fixed by LM)
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000554
Tim Petersd7b5e882001-01-25 03:36:26 +0000555 # Multithreading note: When two threads call this function
556 # simultaneously, it is possible that they will receive the
557 # same return value. The window is very small though. To
558 # avoid this, you have to use a lock around all calls. (I
559 # didn't want to slow this down in the serial case by using a
560 # lock here.)
Guido van Rossumd03e1191998-05-29 17:51:31 +0000561
Tim Petersd7b5e882001-01-25 03:36:26 +0000562 random = self.random
563 z = self.gauss_next
564 self.gauss_next = None
565 if z is None:
566 x2pi = random() * TWOPI
567 g2rad = _sqrt(-2.0 * _log(1.0 - random()))
568 z = _cos(x2pi) * g2rad
569 self.gauss_next = _sin(x2pi) * g2rad
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000570
Tim Petersd7b5e882001-01-25 03:36:26 +0000571 return mu + z*sigma
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000572
Tim Peterscd804102001-01-25 20:25:57 +0000573## -------------------- beta --------------------
Tim Peters85e2e472001-01-26 06:49:56 +0000574## See
575## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470
576## for Ivan Frohne's insightful analysis of why the original implementation:
577##
578## def betavariate(self, alpha, beta):
579## # Discrete Event Simulation in C, pp 87-88.
580##
581## y = self.expovariate(alpha)
582## z = self.expovariate(1.0/beta)
583## return z/(y+z)
584##
585## was dead wrong, and how it probably got that way.
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000586
Tim Petersd7b5e882001-01-25 03:36:26 +0000587 def betavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000588 """Beta distribution.
589
Thomas Woutersb2137042007-02-01 18:02:27 +0000590 Conditions on the parameters are alpha > 0 and beta > 0.
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000591 Returned values range between 0 and 1.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000592
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000593 """
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000594
Tim Peters85e2e472001-01-26 06:49:56 +0000595 # This version due to Janne Sinkkonen, and matches all the std
596 # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
597 y = self.gammavariate(alpha, 1.)
598 if y == 0:
599 return 0.0
600 else:
601 return y / (y + self.gammavariate(beta, 1.))
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000602
Tim Peterscd804102001-01-25 20:25:57 +0000603## -------------------- Pareto --------------------
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000604
Tim Petersd7b5e882001-01-25 03:36:26 +0000605 def paretovariate(self, alpha):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000606 """Pareto distribution. alpha is the shape parameter."""
Tim Petersd7b5e882001-01-25 03:36:26 +0000607 # Jain, pg. 495
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000608
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000609 u = 1.0 - self.random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000610 return 1.0 / pow(u, 1.0/alpha)
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000611
Tim Peterscd804102001-01-25 20:25:57 +0000612## -------------------- Weibull --------------------
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000613
Tim Petersd7b5e882001-01-25 03:36:26 +0000614 def weibullvariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000615 """Weibull distribution.
616
617 alpha is the scale parameter and beta is the shape parameter.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000618
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000619 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000620 # Jain, pg. 499; bug fix courtesy Bill Arms
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000621
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000622 u = 1.0 - self.random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000623 return alpha * pow(-_log(u), 1.0/beta)
Guido van Rossum6c395ba1999-08-18 13:53:28 +0000624
Raymond Hettinger23f12412004-09-13 22:23:21 +0000625## --------------- Operating System Random Source ------------------
Raymond Hettinger356a4592004-08-30 06:14:31 +0000626
Raymond Hettinger23f12412004-09-13 22:23:21 +0000627class SystemRandom(Random):
628 """Alternate random number generator using sources provided
629 by the operating system (such as /dev/urandom on Unix or
630 CryptGenRandom on Windows).
Raymond Hettinger356a4592004-08-30 06:14:31 +0000631
632 Not available on all systems (see os.urandom() for details).
633 """
634
635 def random(self):
636 """Get the next random number in the range [0.0, 1.0)."""
Guido van Rossume2a383d2007-01-15 16:59:06 +0000637 return (int(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF
Raymond Hettinger356a4592004-08-30 06:14:31 +0000638
639 def getrandbits(self, k):
640 """getrandbits(k) -> x. Generates a long int with k random bits."""
Raymond Hettinger356a4592004-08-30 06:14:31 +0000641 if k <= 0:
642 raise ValueError('number of bits must be greater than zero')
643 if k != int(k):
644 raise TypeError('number of bits should be an integer')
645 bytes = (k + 7) // 8 # bits / 8 and rounded up
Guido van Rossume2a383d2007-01-15 16:59:06 +0000646 x = int(_hexlify(_urandom(bytes)), 16)
Raymond Hettinger356a4592004-08-30 06:14:31 +0000647 return x >> (bytes * 8 - k) # trim excess bits
648
Raymond Hettinger28de64f2008-01-13 23:40:30 +0000649 def seed(self, *args, **kwds):
Raymond Hettinger23f12412004-09-13 22:23:21 +0000650 "Stub method. Not used for a system random number generator."
Raymond Hettinger356a4592004-08-30 06:14:31 +0000651 return None
Raymond Hettinger356a4592004-08-30 06:14:31 +0000652
653 def _notimplemented(self, *args, **kwds):
Raymond Hettinger23f12412004-09-13 22:23:21 +0000654 "Method should not be called for a system random number generator."
655 raise NotImplementedError('System entropy source does not have state.')
Raymond Hettinger356a4592004-08-30 06:14:31 +0000656 getstate = setstate = _notimplemented
657
Tim Peterscd804102001-01-25 20:25:57 +0000658## -------------------- test program --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000659
Raymond Hettinger62297132003-08-30 01:24:19 +0000660def _test_generator(n, func, args):
Tim Peters0c9886d2001-01-15 01:18:21 +0000661 import time
Guido van Rossumbe19ed72007-02-09 05:37:30 +0000662 print(n, 'times', func.__name__)
Raymond Hettingerb98154e2003-05-24 17:26:02 +0000663 total = 0.0
Tim Peters0c9886d2001-01-15 01:18:21 +0000664 sqsum = 0.0
665 smallest = 1e10
666 largest = -1e10
667 t0 = time.time()
668 for i in range(n):
Raymond Hettinger62297132003-08-30 01:24:19 +0000669 x = func(*args)
Raymond Hettingerb98154e2003-05-24 17:26:02 +0000670 total += x
Tim Peters0c9886d2001-01-15 01:18:21 +0000671 sqsum = sqsum + x*x
672 smallest = min(x, smallest)
673 largest = max(x, largest)
674 t1 = time.time()
Guido van Rossumbe19ed72007-02-09 05:37:30 +0000675 print(round(t1-t0, 3), 'sec,', end=' ')
Raymond Hettingerb98154e2003-05-24 17:26:02 +0000676 avg = total/n
Tim Petersd7b5e882001-01-25 03:36:26 +0000677 stddev = _sqrt(sqsum/n - avg*avg)
Guido van Rossumbe19ed72007-02-09 05:37:30 +0000678 print('avg %g, stddev %g, min %g, max %g' % \
679 (avg, stddev, smallest, largest))
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000680
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000681
682def _test(N=2000):
Raymond Hettinger62297132003-08-30 01:24:19 +0000683 _test_generator(N, random, ())
684 _test_generator(N, normalvariate, (0.0, 1.0))
685 _test_generator(N, lognormvariate, (0.0, 1.0))
686 _test_generator(N, vonmisesvariate, (0.0, 1.0))
687 _test_generator(N, gammavariate, (0.01, 1.0))
688 _test_generator(N, gammavariate, (0.1, 1.0))
689 _test_generator(N, gammavariate, (0.1, 2.0))
690 _test_generator(N, gammavariate, (0.5, 1.0))
691 _test_generator(N, gammavariate, (0.9, 1.0))
692 _test_generator(N, gammavariate, (1.0, 1.0))
693 _test_generator(N, gammavariate, (2.0, 1.0))
694 _test_generator(N, gammavariate, (20.0, 1.0))
695 _test_generator(N, gammavariate, (200.0, 1.0))
696 _test_generator(N, gauss, (0.0, 1.0))
697 _test_generator(N, betavariate, (3.0, 3.0))
Christian Heimesfe337bf2008-03-23 21:54:12 +0000698 _test_generator(N, triangular, (0.0, 1.0, 1.0/3.0))
Tim Peterscd804102001-01-25 20:25:57 +0000699
Tim Peters715c4c42001-01-26 22:56:56 +0000700# Create one instance, seeded from current time, and export its methods
Raymond Hettinger40f62172002-12-29 23:03:38 +0000701# as module-level functions. The functions share state across all uses
702#(both in the user's code and in the Python libraries), but that's fine
703# for most programs and is easier for the casual user than making them
704# instantiate their own Random() instance.
705
Tim Petersd7b5e882001-01-25 03:36:26 +0000706_inst = Random()
707seed = _inst.seed
708random = _inst.random
709uniform = _inst.uniform
Christian Heimesfe337bf2008-03-23 21:54:12 +0000710triangular = _inst.triangular
Tim Petersd7b5e882001-01-25 03:36:26 +0000711randint = _inst.randint
712choice = _inst.choice
713randrange = _inst.randrange
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000714sample = _inst.sample
Tim Petersd7b5e882001-01-25 03:36:26 +0000715shuffle = _inst.shuffle
716normalvariate = _inst.normalvariate
717lognormvariate = _inst.lognormvariate
Tim Petersd7b5e882001-01-25 03:36:26 +0000718expovariate = _inst.expovariate
719vonmisesvariate = _inst.vonmisesvariate
720gammavariate = _inst.gammavariate
Tim Petersd7b5e882001-01-25 03:36:26 +0000721gauss = _inst.gauss
722betavariate = _inst.betavariate
723paretovariate = _inst.paretovariate
724weibullvariate = _inst.weibullvariate
725getstate = _inst.getstate
726setstate = _inst.setstate
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000727getrandbits = _inst.getrandbits
Tim Petersd7b5e882001-01-25 03:36:26 +0000728
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000729if __name__ == '__main__':
Tim Petersd7b5e882001-01-25 03:36:26 +0000730 _test()