blob: 568de88e636dacae7dc6815f0c624b2a5d27a8ac [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
Raymond Hettingerbbc50ea2008-03-23 13:32:32 +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.
Tim Peters0e115952006-06-10 22:51:45 +000033* It is one of the most extensively tested generators in existence.
34* Without a direct way to compute N steps forward, the semantics of
35 jumpahead(n) are weakened to simply jump to another distant state and rely
36 on the large period to avoid overlapping sequences.
37* The random() method is implemented in C, executes in a single Python step,
38 and is, therefore, threadsafe.
Tim Peterse360d952001-01-26 10:00:39 +000039
Guido van Rossume7b146f2000-02-04 15:28:42 +000040"""
Guido van Rossumd03e1191998-05-29 17:51:31 +000041
Raymond Hettingerc4f7bab2008-03-23 19:37:53 +000042from __future__ import division
Raymond Hettinger2f726e92003-10-05 09:09:15 +000043from warnings import warn as _warn
44from types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType
Raymond Hettinger91e27c22005-08-19 01:36:35 +000045from math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil
Tim Petersd7b5e882001-01-25 03:36:26 +000046from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
Raymond Hettingerc1c43ca2004-09-05 00:00:42 +000047from os import urandom as _urandom
48from binascii import hexlify as _hexlify
Guido van Rossumff03b1a1994-03-09 12:55:02 +000049
Raymond Hettingerf24eb352002-11-12 17:41:57 +000050__all__ = ["Random","seed","random","uniform","randint","choice","sample",
Skip Montanaro0de65802001-02-15 22:15:14 +000051 "randrange","shuffle","normalvariate","lognormvariate",
Raymond Hettingerbbc50ea2008-03-23 13:32:32 +000052 "expovariate","vonmisesvariate","gammavariate","triangular",
Raymond Hettingerf8a52d32003-08-05 12:23:19 +000053 "gauss","betavariate","paretovariate","weibullvariate",
Raymond Hettinger356a4592004-08-30 06:14:31 +000054 "getstate","setstate","jumpahead", "WichmannHill", "getrandbits",
Raymond Hettinger23f12412004-09-13 22:23:21 +000055 "SystemRandom"]
Tim Petersd7b5e882001-01-25 03:36:26 +000056
57NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
Tim Petersd7b5e882001-01-25 03:36:26 +000058TWOPI = 2.0*_pi
Tim Petersd7b5e882001-01-25 03:36:26 +000059LOG4 = _log(4.0)
Tim Petersd7b5e882001-01-25 03:36:26 +000060SG_MAGICCONST = 1.0 + _log(4.5)
Raymond Hettinger2f726e92003-10-05 09:09:15 +000061BPF = 53 # Number of bits in a float
Tim Peters7c2a85b2004-08-31 02:19:55 +000062RECIP_BPF = 2**-BPF
Tim Petersd7b5e882001-01-25 03:36:26 +000063
Raymond Hettinger356a4592004-08-30 06:14:31 +000064
Tim Petersd7b5e882001-01-25 03:36:26 +000065# Translated by Guido van Rossum from C source provided by
Raymond Hettinger40f62172002-12-29 23:03:38 +000066# Adrian Baddeley. Adapted by Raymond Hettinger for use with
Raymond Hettinger3fa19d72004-08-31 01:05:15 +000067# the Mersenne Twister and os.urandom() core generators.
Tim Petersd7b5e882001-01-25 03:36:26 +000068
Raymond Hettinger145a4a02003-01-07 10:25:55 +000069import _random
Raymond Hettinger40f62172002-12-29 23:03:38 +000070
Raymond Hettinger145a4a02003-01-07 10:25:55 +000071class Random(_random.Random):
Raymond Hettingerc32f0332002-05-23 19:44:49 +000072 """Random number generator base class used by bound module functions.
73
74 Used to instantiate instances of Random to get generators that don't
75 share state. Especially useful for multi-threaded programs, creating
76 a different instance of Random for each thread, and using the jumpahead()
77 method to ensure that the generated sequences seen by each thread don't
78 overlap.
79
80 Class Random can also be subclassed if you want to use a different basic
81 generator of your own devising: in that case, override the following
Benjamin Petersonf2eb2b42008-07-30 13:46:53 +000082 methods: random(), seed(), getstate(), setstate() and jumpahead().
83 Optionally, implement a getrandbits() method so that randrange() can cover
84 arbitrarily large ranges.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +000085
Raymond Hettingerc32f0332002-05-23 19:44:49 +000086 """
Tim Petersd7b5e882001-01-25 03:36:26 +000087
Martin v. Löwis6b449f42007-12-03 19:20:02 +000088 VERSION = 3 # used by getstate/setstate
Tim Petersd7b5e882001-01-25 03:36:26 +000089
90 def __init__(self, x=None):
91 """Initialize an instance.
92
93 Optional argument x controls seeding, as for Random.seed().
94 """
95
96 self.seed(x)
Raymond Hettinger40f62172002-12-29 23:03:38 +000097 self.gauss_next = None
Tim Petersd7b5e882001-01-25 03:36:26 +000098
Tim Peters0de88fc2001-02-01 04:59:18 +000099 def seed(self, a=None):
100 """Initialize internal state from hashable object.
Tim Petersd7b5e882001-01-25 03:36:26 +0000101
Raymond Hettinger23f12412004-09-13 22:23:21 +0000102 None or no argument seeds from current time or from an operating
103 system specific randomness source if available.
Tim Peters0de88fc2001-02-01 04:59:18 +0000104
Tim Petersbcd725f2001-02-01 10:06:53 +0000105 If a is not None or an int or long, hash(a) is used instead.
Tim Petersd7b5e882001-01-25 03:36:26 +0000106 """
107
Raymond Hettinger3081d592003-08-09 18:30:57 +0000108 if a is None:
Raymond Hettingerc1c43ca2004-09-05 00:00:42 +0000109 try:
110 a = long(_hexlify(_urandom(16)), 16)
111 except NotImplementedError:
Raymond Hettinger356a4592004-08-30 06:14:31 +0000112 import time
113 a = long(time.time() * 256) # use fractional seconds
Raymond Hettinger356a4592004-08-30 06:14:31 +0000114
Raymond Hettinger145a4a02003-01-07 10:25:55 +0000115 super(Random, self).seed(a)
Tim Peters46c04e12002-05-05 20:40:00 +0000116 self.gauss_next = None
117
Tim Peterscd804102001-01-25 20:25:57 +0000118 def getstate(self):
119 """Return internal state; can be passed to setstate() later."""
Raymond Hettinger145a4a02003-01-07 10:25:55 +0000120 return self.VERSION, super(Random, self).getstate(), self.gauss_next
Tim Peterscd804102001-01-25 20:25:57 +0000121
122 def setstate(self, state):
123 """Restore internal state from object returned by getstate()."""
124 version = state[0]
Martin v. Löwis6b449f42007-12-03 19:20:02 +0000125 if version == 3:
Raymond Hettinger40f62172002-12-29 23:03:38 +0000126 version, internalstate, self.gauss_next = state
Raymond Hettinger145a4a02003-01-07 10:25:55 +0000127 super(Random, self).setstate(internalstate)
Martin v. Löwis6b449f42007-12-03 19:20:02 +0000128 elif version == 2:
129 version, internalstate, self.gauss_next = state
130 # In version 2, the state was saved as signed ints, which causes
131 # inconsistencies between 32/64-bit systems. The state is
132 # really unsigned 32-bit ints, so we convert negative ints from
133 # version 2 to positive longs for version 3.
134 try:
135 internalstate = tuple( long(x) % (2**32) for x in internalstate )
136 except ValueError, e:
137 raise TypeError, e
138 super(Random, self).setstate(internalstate)
Tim Peterscd804102001-01-25 20:25:57 +0000139 else:
140 raise ValueError("state with version %s passed to "
141 "Random.setstate() of version %s" %
142 (version, self.VERSION))
143
Tim Peterscd804102001-01-25 20:25:57 +0000144## ---- Methods below this point do not need to be overridden when
145## ---- subclassing for the purpose of using a different core generator.
146
147## -------------------- pickle support -------------------
148
149 def __getstate__(self): # for pickle
150 return self.getstate()
151
152 def __setstate__(self, state): # for pickle
153 self.setstate(state)
154
Raymond Hettinger5f078ff2003-06-24 20:29:04 +0000155 def __reduce__(self):
156 return self.__class__, (), self.getstate()
157
Tim Peterscd804102001-01-25 20:25:57 +0000158## -------------------- integer methods -------------------
159
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000160 def randrange(self, start, stop=None, step=1, int=int, default=None,
161 maxwidth=1L<<BPF):
Tim Petersd7b5e882001-01-25 03:36:26 +0000162 """Choose a random item from range(start, stop[, step]).
163
164 This fixes the problem with randint() which includes the
165 endpoint; in Python this is usually not what you want.
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000166 Do not supply the 'int', 'default', and 'maxwidth' arguments.
Tim Petersd7b5e882001-01-25 03:36:26 +0000167 """
168
169 # This code is a bit messy to make it fast for the
Tim Peters9146f272002-08-16 03:41:39 +0000170 # common case while still doing adequate error checking.
Tim Petersd7b5e882001-01-25 03:36:26 +0000171 istart = int(start)
172 if istart != start:
173 raise ValueError, "non-integer arg 1 for randrange()"
174 if stop is default:
175 if istart > 0:
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000176 if istart >= maxwidth:
177 return self._randbelow(istart)
Tim Petersd7b5e882001-01-25 03:36:26 +0000178 return int(self.random() * istart)
179 raise ValueError, "empty range for randrange()"
Tim Peters9146f272002-08-16 03:41:39 +0000180
181 # stop argument supplied.
Tim Petersd7b5e882001-01-25 03:36:26 +0000182 istop = int(stop)
183 if istop != stop:
184 raise ValueError, "non-integer stop for randrange()"
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000185 width = istop - istart
186 if step == 1 and width > 0:
Tim Peters76ca1d42003-06-19 03:46:46 +0000187 # Note that
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000188 # int(istart + self.random()*width)
Tim Peters76ca1d42003-06-19 03:46:46 +0000189 # instead would be incorrect. For example, consider istart
190 # = -2 and istop = 0. Then the guts would be in
191 # -2.0 to 0.0 exclusive on both ends (ignoring that random()
192 # might return 0.0), and because int() truncates toward 0, the
193 # final result would be -1 or 0 (instead of -2 or -1).
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000194 # istart + int(self.random()*width)
Tim Peters76ca1d42003-06-19 03:46:46 +0000195 # would also be incorrect, for a subtler reason: the RHS
196 # can return a long, and then randrange() would also return
197 # a long, but we're supposed to return an int (for backward
198 # compatibility).
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000199
200 if width >= maxwidth:
Tim Peters58eb11c2004-01-18 20:29:55 +0000201 return int(istart + self._randbelow(width))
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000202 return int(istart + int(self.random()*width))
Tim Petersd7b5e882001-01-25 03:36:26 +0000203 if step == 1:
Raymond Hettinger2f726e92003-10-05 09:09:15 +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:
209 raise ValueError, "non-integer step for randrange()"
210 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:
215 raise ValueError, "zero step for randrange()"
216
217 if n <= 0:
218 raise ValueError, "empty range for randrange()"
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000219
220 if n >= maxwidth:
Raymond Hettinger94547f72006-12-20 06:42:06 +0000221 return istart + istep*self._randbelow(n)
Tim Petersd7b5e882001-01-25 03:36:26 +0000222 return istart + istep*int(self.random() * n)
223
224 def randint(self, a, b):
Tim Peterscd804102001-01-25 20:25:57 +0000225 """Return random integer in range [a, b], including both end points.
Tim Petersd7b5e882001-01-25 03:36:26 +0000226 """
227
228 return self.randrange(a, b+1)
229
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000230 def _randbelow(self, n, _log=_log, int=int, _maxwidth=1L<<BPF,
231 _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType):
232 """Return a random int in the range [0,n)
233
234 Handles the case where n has more bits than returned
235 by a single call to the underlying generator.
236 """
237
238 try:
239 getrandbits = self.getrandbits
240 except AttributeError:
241 pass
242 else:
243 # Only call self.getrandbits if the original random() builtin method
244 # has not been overridden or if a new getrandbits() was supplied.
245 # This assures that the two methods correspond.
246 if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method:
247 k = int(1.00001 + _log(n-1, 2.0)) # 2**k > n-1 > 2**(k-2)
248 r = getrandbits(k)
249 while r >= n:
250 r = getrandbits(k)
251 return r
252 if n >= _maxwidth:
253 _warn("Underlying random() generator does not supply \n"
254 "enough bits to choose from a population range this large")
255 return int(self.random() * n)
256
Tim Peterscd804102001-01-25 20:25:57 +0000257## -------------------- sequence methods -------------------
258
Tim Petersd7b5e882001-01-25 03:36:26 +0000259 def choice(self, seq):
260 """Choose a random element from a non-empty sequence."""
Raymond Hettinger5dae5052004-06-07 02:07:15 +0000261 return seq[int(self.random() * len(seq))] # raises IndexError if seq is empty
Tim Petersd7b5e882001-01-25 03:36:26 +0000262
263 def shuffle(self, x, random=None, int=int):
264 """x, random=random.random -> shuffle list x in place; return None.
265
266 Optional arg random is a 0-argument function returning a random
267 float in [0.0, 1.0); by default, the standard random.random.
Tim Petersd7b5e882001-01-25 03:36:26 +0000268 """
269
270 if random is None:
271 random = self.random
Raymond Hettinger85c20a42003-11-06 14:06:48 +0000272 for i in reversed(xrange(1, len(x))):
Tim Peterscd804102001-01-25 20:25:57 +0000273 # pick an element in x[:i+1] with which to exchange x[i]
Tim Petersd7b5e882001-01-25 03:36:26 +0000274 j = int(random() * (i+1))
275 x[i], x[j] = x[j], x[i]
276
Raymond Hettingerfdbe5222003-06-13 07:01:51 +0000277 def sample(self, population, k):
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000278 """Chooses k unique random elements from a population sequence.
279
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000280 Returns a new list containing elements from the population while
281 leaving the original population unchanged. The resulting list is
282 in selection order so that all sub-slices will also be valid random
283 samples. This allows raffle winners (the sample) to be partitioned
284 into grand prize and second place winners (the subslices).
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000285
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000286 Members of the population need not be hashable or unique. If the
287 population contains repeats, then each occurrence is a possible
288 selection in the sample.
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000289
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000290 To choose a sample in a range of integers, use xrange as an argument.
291 This is especially fast and space efficient for sampling from a
292 large population: sample(xrange(10000000), 60)
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000293 """
294
Tim Petersc17976e2006-04-01 00:26:53 +0000295 # XXX Although the documentation says `population` is "a sequence",
296 # XXX attempts are made to cater to any iterable with a __len__
297 # XXX method. This has had mixed success. Examples from both
298 # XXX sides: sets work fine, and should become officially supported;
299 # XXX dicts are much harder, and have failed in various subtle
300 # XXX ways across attempts. Support for mapping types should probably
301 # XXX be dropped (and users should pass mapping.keys() or .values()
302 # XXX explicitly).
303
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000304 # Sampling without replacement entails tracking either potential
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000305 # selections (the pool) in a list or previous selections in a set.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000306
Jeremy Hylton2b55d352004-02-23 17:27:57 +0000307 # When the number of selections is small compared to the
308 # population, then tracking selections is efficient, requiring
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000309 # only a small set and an occasional reselection. For
Jeremy Hylton2b55d352004-02-23 17:27:57 +0000310 # a larger number of selections, the pool tracking method is
311 # preferred since the list takes less space than the
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000312 # set and it doesn't suffer from frequent reselections.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000313
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000314 n = len(population)
315 if not 0 <= k <= n:
316 raise ValueError, "sample larger than population"
Raymond Hettinger8b9aa8d2003-01-04 05:20:33 +0000317 random = self.random
Raymond Hettingerfdbe5222003-06-13 07:01:51 +0000318 _int = int
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000319 result = [None] * k
Raymond Hettinger91e27c22005-08-19 01:36:35 +0000320 setsize = 21 # size of a small set minus size of an empty list
321 if k > 5:
Tim Peters9e34c042005-08-26 15:20:46 +0000322 setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets
Tim Petersc17976e2006-04-01 00:26:53 +0000323 if n <= setsize or hasattr(population, "keys"):
324 # An n-length list is smaller than a k-length set, or this is a
325 # mapping type so the other algorithm wouldn't work.
Raymond Hettinger311f4192002-11-18 09:01:24 +0000326 pool = list(population)
327 for i in xrange(k): # invariant: non-selected at [0,n-i)
Raymond Hettingerfdbe5222003-06-13 07:01:51 +0000328 j = _int(random() * (n-i))
Raymond Hettinger311f4192002-11-18 09:01:24 +0000329 result[i] = pool[j]
Raymond Hettinger8b9aa8d2003-01-04 05:20:33 +0000330 pool[j] = pool[n-i-1] # move non-selected item into vacancy
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000331 else:
Raymond Hettinger66d09f12003-09-06 04:25:54 +0000332 try:
Raymond Hettinger3c3346d2006-03-29 09:13:13 +0000333 selected = set()
334 selected_add = selected.add
335 for i in xrange(k):
Raymond Hettingerfdbe5222003-06-13 07:01:51 +0000336 j = _int(random() * n)
Raymond Hettinger3c3346d2006-03-29 09:13:13 +0000337 while j in selected:
338 j = _int(random() * n)
339 selected_add(j)
340 result[i] = population[j]
Tim Petersc17976e2006-04-01 00:26:53 +0000341 except (TypeError, KeyError): # handle (at least) sets
Raymond Hettinger3c3346d2006-03-29 09:13:13 +0000342 if isinstance(population, list):
343 raise
Tim Petersc17976e2006-04-01 00:26:53 +0000344 return self.sample(tuple(population), k)
Raymond Hettinger311f4192002-11-18 09:01:24 +0000345 return result
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000346
Tim Peterscd804102001-01-25 20:25:57 +0000347## -------------------- real-valued distributions -------------------
348
349## -------------------- uniform distribution -------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000350
351 def uniform(self, a, b):
Raymond Hettinger2c0cdca2009-06-11 23:14:53 +0000352 "Get a random number in the range [a, b) or [a, b] depending on rounding."
Tim Petersd7b5e882001-01-25 03:36:26 +0000353 return a + (b-a) * self.random()
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000354
Raymond Hettingerbbc50ea2008-03-23 13:32:32 +0000355## -------------------- triangular --------------------
356
Raymond Hettingerc4f7bab2008-03-23 19:37:53 +0000357 def triangular(self, low=0.0, high=1.0, mode=None):
Raymond Hettingerbbc50ea2008-03-23 13:32:32 +0000358 """Triangular distribution.
359
360 Continuous distribution bounded by given lower and upper limits,
361 and having a given mode value in-between.
362
363 http://en.wikipedia.org/wiki/Triangular_distribution
364
365 """
366 u = self.random()
Raymond Hettingerc4f7bab2008-03-23 19:37:53 +0000367 c = 0.5 if mode is None else (mode - low) / (high - low)
Raymond Hettingerbbc50ea2008-03-23 13:32:32 +0000368 if u > c:
Raymond Hettingerc4f7bab2008-03-23 19:37:53 +0000369 u = 1.0 - u
370 c = 1.0 - c
Raymond Hettingerbbc50ea2008-03-23 13:32:32 +0000371 low, high = high, low
372 return low + (high - low) * (u * c) ** 0.5
373
Tim Peterscd804102001-01-25 20:25:57 +0000374## -------------------- normal distribution --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000375
Tim Petersd7b5e882001-01-25 03:36:26 +0000376 def normalvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000377 """Normal distribution.
378
379 mu is the mean, and sigma is the standard deviation.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000380
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000381 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000382 # mu = mean, sigma = standard deviation
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000383
Tim Petersd7b5e882001-01-25 03:36:26 +0000384 # Uses Kinderman and Monahan method. Reference: Kinderman,
385 # A.J. and Monahan, J.F., "Computer generation of random
386 # variables using the ratio of uniform deviates", ACM Trans
387 # Math Software, 3, (1977), pp257-260.
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000388
Tim Petersd7b5e882001-01-25 03:36:26 +0000389 random = self.random
Raymond Hettinger42406e62005-04-30 09:02:51 +0000390 while 1:
Tim Peters0c9886d2001-01-15 01:18:21 +0000391 u1 = random()
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000392 u2 = 1.0 - random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000393 z = NV_MAGICCONST*(u1-0.5)/u2
394 zz = z*z/4.0
395 if zz <= -_log(u2):
396 break
397 return mu + z*sigma
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000398
Tim Peterscd804102001-01-25 20:25:57 +0000399## -------------------- lognormal distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000400
401 def lognormvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000402 """Log normal distribution.
403
404 If you take the natural logarithm of this distribution, you'll get a
405 normal distribution with mean mu and standard deviation sigma.
406 mu can have any value, and sigma must be greater than zero.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000407
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000408 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000409 return _exp(self.normalvariate(mu, sigma))
410
Tim Peterscd804102001-01-25 20:25:57 +0000411## -------------------- exponential distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000412
413 def expovariate(self, lambd):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000414 """Exponential distribution.
415
Mark Dickinsone6dc5312009-01-07 17:48:33 +0000416 lambd is 1.0 divided by the desired mean. It should be
417 nonzero. (The parameter would be called "lambda", but that is
418 a reserved word in Python.) Returned values range from 0 to
419 positive infinity if lambd is positive, and from negative
420 infinity to 0 if lambd is negative.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000421
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000422 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000423 # lambd: rate lambd = 1/mean
424 # ('lambda' is a Python reserved word)
425
426 random = self.random
Tim Peters0c9886d2001-01-15 01:18:21 +0000427 u = random()
428 while u <= 1e-7:
429 u = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000430 return -_log(u)/lambd
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000431
Tim Peterscd804102001-01-25 20:25:57 +0000432## -------------------- von Mises distribution --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000433
Tim Petersd7b5e882001-01-25 03:36:26 +0000434 def vonmisesvariate(self, mu, kappa):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000435 """Circular data distribution.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000436
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000437 mu is the mean angle, expressed in radians between 0 and 2*pi, and
438 kappa is the concentration parameter, which must be greater than or
439 equal to zero. If kappa is equal to zero, this distribution reduces
440 to a uniform random angle over the range 0 to 2*pi.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000441
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000442 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000443 # mu: mean angle (in radians between 0 and 2*pi)
444 # kappa: concentration parameter kappa (>= 0)
445 # if kappa = 0 generate uniform random angle
446
447 # Based upon an algorithm published in: Fisher, N.I.,
448 # "Statistical Analysis of Circular Data", Cambridge
449 # University Press, 1993.
450
451 # Thanks to Magnus Kessler for a correction to the
452 # implementation of step 4.
453
454 random = self.random
455 if kappa <= 1e-6:
456 return TWOPI * random()
457
458 a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa)
459 b = (a - _sqrt(2.0 * a))/(2.0 * kappa)
460 r = (1.0 + b * b)/(2.0 * b)
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000461
Raymond Hettinger42406e62005-04-30 09:02:51 +0000462 while 1:
Tim Peters0c9886d2001-01-15 01:18:21 +0000463 u1 = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000464
465 z = _cos(_pi * u1)
466 f = (1.0 + r * z)/(r + z)
467 c = kappa * (r - f)
468
469 u2 = random()
470
Raymond Hettinger42406e62005-04-30 09:02:51 +0000471 if u2 < c * (2.0 - c) or u2 <= c * _exp(1.0 - c):
Tim Peters0c9886d2001-01-15 01:18:21 +0000472 break
Tim Petersd7b5e882001-01-25 03:36:26 +0000473
474 u3 = random()
475 if u3 > 0.5:
476 theta = (mu % TWOPI) + _acos(f)
477 else:
478 theta = (mu % TWOPI) - _acos(f)
479
480 return theta
481
Tim Peterscd804102001-01-25 20:25:57 +0000482## -------------------- gamma distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000483
484 def gammavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000485 """Gamma distribution. Not the gamma function!
486
487 Conditions on the parameters are alpha > 0 and beta > 0.
488
489 """
Tim Peters8ac14952002-05-23 15:15:30 +0000490
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000491 # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
Tim Peters8ac14952002-05-23 15:15:30 +0000492
Guido van Rossum570764d2002-05-14 14:08:12 +0000493 # Warning: a few older sources define the gamma distribution in terms
494 # of alpha > -1.0
495 if alpha <= 0.0 or beta <= 0.0:
496 raise ValueError, 'gammavariate: alpha and beta must be > 0.0'
Tim Peters8ac14952002-05-23 15:15:30 +0000497
Tim Petersd7b5e882001-01-25 03:36:26 +0000498 random = self.random
Tim Petersd7b5e882001-01-25 03:36:26 +0000499 if alpha > 1.0:
500
501 # Uses R.C.H. Cheng, "The generation of Gamma
502 # variables with non-integral shape parameters",
503 # Applied Statistics, (1977), 26, No. 1, p71-74
504
Raymond Hettingerca6cdc22002-05-13 23:40:14 +0000505 ainv = _sqrt(2.0 * alpha - 1.0)
506 bbb = alpha - LOG4
507 ccc = alpha + ainv
Tim Peters8ac14952002-05-23 15:15:30 +0000508
Raymond Hettinger42406e62005-04-30 09:02:51 +0000509 while 1:
Tim Petersd7b5e882001-01-25 03:36:26 +0000510 u1 = random()
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000511 if not 1e-7 < u1 < .9999999:
512 continue
513 u2 = 1.0 - random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000514 v = _log(u1/(1.0-u1))/ainv
515 x = alpha*_exp(v)
516 z = u1*u1*u2
517 r = bbb+ccc*v-x
518 if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z):
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000519 return x * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000520
521 elif alpha == 1.0:
522 # expovariate(1)
523 u = random()
524 while u <= 1e-7:
525 u = random()
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000526 return -_log(u) * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000527
528 else: # alpha is between 0 and 1 (exclusive)
529
530 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
531
Raymond Hettinger42406e62005-04-30 09:02:51 +0000532 while 1:
Tim Petersd7b5e882001-01-25 03:36:26 +0000533 u = random()
534 b = (_e + alpha)/_e
535 p = b*u
536 if p <= 1.0:
Raymond Hettinger42406e62005-04-30 09:02:51 +0000537 x = p ** (1.0/alpha)
Tim Petersd7b5e882001-01-25 03:36:26 +0000538 else:
Tim Petersd7b5e882001-01-25 03:36:26 +0000539 x = -_log((b-p)/alpha)
540 u1 = random()
Raymond Hettinger42406e62005-04-30 09:02:51 +0000541 if p > 1.0:
542 if u1 <= x ** (alpha - 1.0):
543 break
544 elif u1 <= _exp(-x):
Tim Petersd7b5e882001-01-25 03:36:26 +0000545 break
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000546 return x * beta
547
Tim Peterscd804102001-01-25 20:25:57 +0000548## -------------------- Gauss (faster alternative) --------------------
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000549
Tim Petersd7b5e882001-01-25 03:36:26 +0000550 def gauss(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000551 """Gaussian distribution.
552
553 mu is the mean, and sigma is the standard deviation. This is
554 slightly faster than the normalvariate() function.
555
556 Not thread-safe without a lock around calls.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000557
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000558 """
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000559
Tim Petersd7b5e882001-01-25 03:36:26 +0000560 # When x and y are two variables from [0, 1), uniformly
561 # distributed, then
562 #
563 # cos(2*pi*x)*sqrt(-2*log(1-y))
564 # sin(2*pi*x)*sqrt(-2*log(1-y))
565 #
566 # are two *independent* variables with normal distribution
567 # (mu = 0, sigma = 1).
568 # (Lambert Meertens)
569 # (corrected version; bug discovered by Mike Miller, fixed by LM)
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000570
Tim Petersd7b5e882001-01-25 03:36:26 +0000571 # Multithreading note: When two threads call this function
572 # simultaneously, it is possible that they will receive the
573 # same return value. The window is very small though. To
574 # avoid this, you have to use a lock around all calls. (I
575 # didn't want to slow this down in the serial case by using a
576 # lock here.)
Guido van Rossumd03e1191998-05-29 17:51:31 +0000577
Tim Petersd7b5e882001-01-25 03:36:26 +0000578 random = self.random
579 z = self.gauss_next
580 self.gauss_next = None
581 if z is None:
582 x2pi = random() * TWOPI
583 g2rad = _sqrt(-2.0 * _log(1.0 - random()))
584 z = _cos(x2pi) * g2rad
585 self.gauss_next = _sin(x2pi) * g2rad
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000586
Tim Petersd7b5e882001-01-25 03:36:26 +0000587 return mu + z*sigma
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000588
Tim Peterscd804102001-01-25 20:25:57 +0000589## -------------------- beta --------------------
Tim Peters85e2e472001-01-26 06:49:56 +0000590## See
591## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470
592## for Ivan Frohne's insightful analysis of why the original implementation:
593##
594## def betavariate(self, alpha, beta):
595## # Discrete Event Simulation in C, pp 87-88.
596##
597## y = self.expovariate(alpha)
598## z = self.expovariate(1.0/beta)
599## return z/(y+z)
600##
601## was dead wrong, and how it probably got that way.
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000602
Tim Petersd7b5e882001-01-25 03:36:26 +0000603 def betavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000604 """Beta distribution.
605
Raymond Hettinger1b0ce852007-01-19 18:07:18 +0000606 Conditions on the parameters are alpha > 0 and beta > 0.
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000607 Returned values range between 0 and 1.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000608
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000609 """
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000610
Tim Peters85e2e472001-01-26 06:49:56 +0000611 # This version due to Janne Sinkkonen, and matches all the std
612 # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
613 y = self.gammavariate(alpha, 1.)
614 if y == 0:
615 return 0.0
616 else:
617 return y / (y + self.gammavariate(beta, 1.))
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000618
Tim Peterscd804102001-01-25 20:25:57 +0000619## -------------------- Pareto --------------------
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000620
Tim Petersd7b5e882001-01-25 03:36:26 +0000621 def paretovariate(self, alpha):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000622 """Pareto distribution. alpha is the shape parameter."""
Tim Petersd7b5e882001-01-25 03:36:26 +0000623 # Jain, pg. 495
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000624
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000625 u = 1.0 - self.random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000626 return 1.0 / pow(u, 1.0/alpha)
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000627
Tim Peterscd804102001-01-25 20:25:57 +0000628## -------------------- Weibull --------------------
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000629
Tim Petersd7b5e882001-01-25 03:36:26 +0000630 def weibullvariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000631 """Weibull distribution.
632
633 alpha is the scale parameter and beta is the shape parameter.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000634
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000635 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000636 # Jain, pg. 499; bug fix courtesy Bill Arms
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000637
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000638 u = 1.0 - self.random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000639 return alpha * pow(-_log(u), 1.0/beta)
Guido van Rossum6c395ba1999-08-18 13:53:28 +0000640
Raymond Hettinger40f62172002-12-29 23:03:38 +0000641## -------------------- Wichmann-Hill -------------------
642
643class WichmannHill(Random):
644
645 VERSION = 1 # used by getstate/setstate
646
647 def seed(self, a=None):
648 """Initialize internal state from hashable object.
649
Raymond Hettinger23f12412004-09-13 22:23:21 +0000650 None or no argument seeds from current time or from an operating
651 system specific randomness source if available.
Raymond Hettinger40f62172002-12-29 23:03:38 +0000652
653 If a is not None or an int or long, hash(a) is used instead.
654
655 If a is an int or long, a is used directly. Distinct values between
656 0 and 27814431486575L inclusive are guaranteed to yield distinct
657 internal states (this guarantee is specific to the default
658 Wichmann-Hill generator).
659 """
660
661 if a is None:
Raymond Hettingerc1c43ca2004-09-05 00:00:42 +0000662 try:
663 a = long(_hexlify(_urandom(16)), 16)
664 except NotImplementedError:
Raymond Hettinger356a4592004-08-30 06:14:31 +0000665 import time
666 a = long(time.time() * 256) # use fractional seconds
Raymond Hettinger40f62172002-12-29 23:03:38 +0000667
668 if not isinstance(a, (int, long)):
669 a = hash(a)
670
671 a, x = divmod(a, 30268)
672 a, y = divmod(a, 30306)
673 a, z = divmod(a, 30322)
674 self._seed = int(x)+1, int(y)+1, int(z)+1
675
676 self.gauss_next = None
677
678 def random(self):
679 """Get the next random number in the range [0.0, 1.0)."""
680
681 # Wichman-Hill random number generator.
682 #
683 # Wichmann, B. A. & Hill, I. D. (1982)
684 # Algorithm AS 183:
685 # An efficient and portable pseudo-random number generator
686 # Applied Statistics 31 (1982) 188-190
687 #
688 # see also:
689 # Correction to Algorithm AS 183
690 # Applied Statistics 33 (1984) 123
691 #
692 # McLeod, A. I. (1985)
693 # A remark on Algorithm AS 183
694 # Applied Statistics 34 (1985),198-200
695
696 # This part is thread-unsafe:
697 # BEGIN CRITICAL SECTION
698 x, y, z = self._seed
699 x = (171 * x) % 30269
700 y = (172 * y) % 30307
701 z = (170 * z) % 30323
702 self._seed = x, y, z
703 # END CRITICAL SECTION
704
705 # Note: on a platform using IEEE-754 double arithmetic, this can
706 # never return 0.0 (asserted by Tim; proof too long for a comment).
707 return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0
708
709 def getstate(self):
710 """Return internal state; can be passed to setstate() later."""
711 return self.VERSION, self._seed, self.gauss_next
712
713 def setstate(self, state):
714 """Restore internal state from object returned by getstate()."""
715 version = state[0]
716 if version == 1:
717 version, self._seed, self.gauss_next = state
718 else:
719 raise ValueError("state with version %s passed to "
720 "Random.setstate() of version %s" %
721 (version, self.VERSION))
722
723 def jumpahead(self, n):
724 """Act as if n calls to random() were made, but quickly.
725
726 n is an int, greater than or equal to 0.
727
728 Example use: If you have 2 threads and know that each will
729 consume no more than a million random numbers, create two Random
730 objects r1 and r2, then do
731 r2.setstate(r1.getstate())
732 r2.jumpahead(1000000)
733 Then r1 and r2 will use guaranteed-disjoint segments of the full
734 period.
735 """
736
737 if not n >= 0:
738 raise ValueError("n must be >= 0")
739 x, y, z = self._seed
740 x = int(x * pow(171, n, 30269)) % 30269
741 y = int(y * pow(172, n, 30307)) % 30307
742 z = int(z * pow(170, n, 30323)) % 30323
743 self._seed = x, y, z
744
745 def __whseed(self, x=0, y=0, z=0):
746 """Set the Wichmann-Hill seed from (x, y, z).
747
748 These must be integers in the range [0, 256).
749 """
750
751 if not type(x) == type(y) == type(z) == int:
752 raise TypeError('seeds must be integers')
753 if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256):
754 raise ValueError('seeds must be in range(0, 256)')
755 if 0 == x == y == z:
756 # Initialize from current time
757 import time
758 t = long(time.time() * 256)
759 t = int((t&0xffffff) ^ (t>>24))
760 t, x = divmod(t, 256)
761 t, y = divmod(t, 256)
762 t, z = divmod(t, 256)
763 # Zero is a poor seed, so substitute 1
764 self._seed = (x or 1, y or 1, z or 1)
765
766 self.gauss_next = None
767
768 def whseed(self, a=None):
769 """Seed from hashable object's hash code.
770
771 None or no argument seeds from current time. It is not guaranteed
772 that objects with distinct hash codes lead to distinct internal
773 states.
774
775 This is obsolete, provided for compatibility with the seed routine
776 used prior to Python 2.1. Use the .seed() method instead.
777 """
778
779 if a is None:
780 self.__whseed()
781 return
782 a = hash(a)
783 a, x = divmod(a, 256)
784 a, y = divmod(a, 256)
785 a, z = divmod(a, 256)
786 x = (x + a) % 256 or 1
787 y = (y + a) % 256 or 1
788 z = (z + a) % 256 or 1
789 self.__whseed(x, y, z)
790
Raymond Hettinger23f12412004-09-13 22:23:21 +0000791## --------------- Operating System Random Source ------------------
Raymond Hettinger356a4592004-08-30 06:14:31 +0000792
Raymond Hettinger23f12412004-09-13 22:23:21 +0000793class SystemRandom(Random):
794 """Alternate random number generator using sources provided
795 by the operating system (such as /dev/urandom on Unix or
796 CryptGenRandom on Windows).
Raymond Hettinger356a4592004-08-30 06:14:31 +0000797
798 Not available on all systems (see os.urandom() for details).
799 """
800
801 def random(self):
802 """Get the next random number in the range [0.0, 1.0)."""
Tim Peters7c2a85b2004-08-31 02:19:55 +0000803 return (long(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF
Raymond Hettinger356a4592004-08-30 06:14:31 +0000804
805 def getrandbits(self, k):
806 """getrandbits(k) -> x. Generates a long int with k random bits."""
Raymond Hettinger356a4592004-08-30 06:14:31 +0000807 if k <= 0:
808 raise ValueError('number of bits must be greater than zero')
809 if k != int(k):
810 raise TypeError('number of bits should be an integer')
811 bytes = (k + 7) // 8 # bits / 8 and rounded up
812 x = long(_hexlify(_urandom(bytes)), 16)
813 return x >> (bytes * 8 - k) # trim excess bits
814
815 def _stub(self, *args, **kwds):
Raymond Hettinger23f12412004-09-13 22:23:21 +0000816 "Stub method. Not used for a system random number generator."
Raymond Hettinger356a4592004-08-30 06:14:31 +0000817 return None
818 seed = jumpahead = _stub
819
820 def _notimplemented(self, *args, **kwds):
Raymond Hettinger23f12412004-09-13 22:23:21 +0000821 "Method should not be called for a system random number generator."
822 raise NotImplementedError('System entropy source does not have state.')
Raymond Hettinger356a4592004-08-30 06:14:31 +0000823 getstate = setstate = _notimplemented
824
Tim Peterscd804102001-01-25 20:25:57 +0000825## -------------------- test program --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000826
Raymond Hettinger62297132003-08-30 01:24:19 +0000827def _test_generator(n, func, args):
Tim Peters0c9886d2001-01-15 01:18:21 +0000828 import time
Raymond Hettinger62297132003-08-30 01:24:19 +0000829 print n, 'times', func.__name__
Raymond Hettingerb98154e2003-05-24 17:26:02 +0000830 total = 0.0
Tim Peters0c9886d2001-01-15 01:18:21 +0000831 sqsum = 0.0
832 smallest = 1e10
833 largest = -1e10
834 t0 = time.time()
835 for i in range(n):
Raymond Hettinger62297132003-08-30 01:24:19 +0000836 x = func(*args)
Raymond Hettingerb98154e2003-05-24 17:26:02 +0000837 total += x
Tim Peters0c9886d2001-01-15 01:18:21 +0000838 sqsum = sqsum + x*x
839 smallest = min(x, smallest)
840 largest = max(x, largest)
841 t1 = time.time()
842 print round(t1-t0, 3), 'sec,',
Raymond Hettingerb98154e2003-05-24 17:26:02 +0000843 avg = total/n
Tim Petersd7b5e882001-01-25 03:36:26 +0000844 stddev = _sqrt(sqsum/n - avg*avg)
Tim Peters0c9886d2001-01-15 01:18:21 +0000845 print 'avg %g, stddev %g, min %g, max %g' % \
846 (avg, stddev, smallest, largest)
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000847
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000848
849def _test(N=2000):
Raymond Hettinger62297132003-08-30 01:24:19 +0000850 _test_generator(N, random, ())
851 _test_generator(N, normalvariate, (0.0, 1.0))
852 _test_generator(N, lognormvariate, (0.0, 1.0))
853 _test_generator(N, vonmisesvariate, (0.0, 1.0))
854 _test_generator(N, gammavariate, (0.01, 1.0))
855 _test_generator(N, gammavariate, (0.1, 1.0))
856 _test_generator(N, gammavariate, (0.1, 2.0))
857 _test_generator(N, gammavariate, (0.5, 1.0))
858 _test_generator(N, gammavariate, (0.9, 1.0))
859 _test_generator(N, gammavariate, (1.0, 1.0))
860 _test_generator(N, gammavariate, (2.0, 1.0))
861 _test_generator(N, gammavariate, (20.0, 1.0))
862 _test_generator(N, gammavariate, (200.0, 1.0))
863 _test_generator(N, gauss, (0.0, 1.0))
864 _test_generator(N, betavariate, (3.0, 3.0))
Raymond Hettingerbbc50ea2008-03-23 13:32:32 +0000865 _test_generator(N, triangular, (0.0, 1.0, 1.0/3.0))
Tim Peterscd804102001-01-25 20:25:57 +0000866
Tim Peters715c4c42001-01-26 22:56:56 +0000867# Create one instance, seeded from current time, and export its methods
Raymond Hettinger40f62172002-12-29 23:03:38 +0000868# as module-level functions. The functions share state across all uses
869#(both in the user's code and in the Python libraries), but that's fine
870# for most programs and is easier for the casual user than making them
871# instantiate their own Random() instance.
872
Tim Petersd7b5e882001-01-25 03:36:26 +0000873_inst = Random()
874seed = _inst.seed
875random = _inst.random
876uniform = _inst.uniform
Raymond Hettingerbbc50ea2008-03-23 13:32:32 +0000877triangular = _inst.triangular
Tim Petersd7b5e882001-01-25 03:36:26 +0000878randint = _inst.randint
879choice = _inst.choice
880randrange = _inst.randrange
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000881sample = _inst.sample
Tim Petersd7b5e882001-01-25 03:36:26 +0000882shuffle = _inst.shuffle
883normalvariate = _inst.normalvariate
884lognormvariate = _inst.lognormvariate
Tim Petersd7b5e882001-01-25 03:36:26 +0000885expovariate = _inst.expovariate
886vonmisesvariate = _inst.vonmisesvariate
887gammavariate = _inst.gammavariate
Tim Petersd7b5e882001-01-25 03:36:26 +0000888gauss = _inst.gauss
889betavariate = _inst.betavariate
890paretovariate = _inst.paretovariate
891weibullvariate = _inst.weibullvariate
892getstate = _inst.getstate
893setstate = _inst.setstate
Tim Petersd52269b2001-01-25 06:23:18 +0000894jumpahead = _inst.jumpahead
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000895getrandbits = _inst.getrandbits
Tim Petersd7b5e882001-01-25 03:36:26 +0000896
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000897if __name__ == '__main__':
Tim Petersd7b5e882001-01-25 03:36:26 +0000898 _test()