source: python/trunk/Lib/random.py@ 12

Last change on this file since 12 was 2, checked in by Yuri Dario, 15 years ago

Initial import for vendor code.
  • Property svn:eol-style set to native
File size: 31.2 KB
Line 
1"""Random variable generators.
2
3 integers
4 --------
5 uniform within range
6
7 sequences
8 ---------
9 pick random element
10 pick random sample
11 generate random permutation
12
13 distributions on the real line:
14 ------------------------------
15 uniform
16 triangular
17 normal (Gaussian)
18 lognormal
19 negative exponential
20 gamma
21 beta
22 pareto
23 Weibull
24
25 distributions on the circle (angles 0 to 2pi)
26 ---------------------------------------------
27 circular uniform
28 von Mises
29
30General notes on the underlying Mersenne Twister core generator:
31
32* The period is 2**19937-1.
33* 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.
39
40"""
41
42from __future__ import division
43from warnings import warn as _warn
44from types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType
45from math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil
46from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
47from os import urandom as _urandom
48from binascii import hexlify as _hexlify
49
50__all__ = ["Random","seed","random","uniform","randint","choice","sample",
51 "randrange","shuffle","normalvariate","lognormvariate",
52 "expovariate","vonmisesvariate","gammavariate","triangular",
53 "gauss","betavariate","paretovariate","weibullvariate",
54 "getstate","setstate","jumpahead", "WichmannHill", "getrandbits",
55 "SystemRandom"]
56
57NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
58TWOPI = 2.0*_pi
59LOG4 = _log(4.0)
60SG_MAGICCONST = 1.0 + _log(4.5)
61BPF = 53 # Number of bits in a float
62RECIP_BPF = 2**-BPF
63
64
65# Translated by Guido van Rossum from C source provided by
66# Adrian Baddeley. Adapted by Raymond Hettinger for use with
67# the Mersenne Twister and os.urandom() core generators.
68
69import _random
70
71class Random(_random.Random):
72 """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
82 methods: random(), seed(), getstate(), setstate() and jumpahead().
83 Optionally, implement a getrandbits() method so that randrange() can cover
84 arbitrarily large ranges.
85
86 """
87
88 VERSION = 3 # used by getstate/setstate
89
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)
97 self.gauss_next = None
98
99 def seed(self, a=None):
100 """Initialize internal state from hashable object.
101
102 None or no argument seeds from current time or from an operating
103 system specific randomness source if available.
104
105 If a is not None or an int or long, hash(a) is used instead.
106 """
107
108 if a is None:
109 try:
110 a = long(_hexlify(_urandom(16)), 16)
111 except NotImplementedError:
112 import time
113 a = long(time.time() * 256) # use fractional seconds
114
115 super(Random, self).seed(a)
116 self.gauss_next = None
117
118 def getstate(self):
119 """Return internal state; can be passed to setstate() later."""
120 return self.VERSION, super(Random, self).getstate(), self.gauss_next
121
122 def setstate(self, state):
123 """Restore internal state from object returned by getstate()."""
124 version = state[0]
125 if version == 3:
126 version, internalstate, self.gauss_next = state
127 super(Random, self).setstate(internalstate)
128 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)
139 else:
140 raise ValueError("state with version %s passed to "
141 "Random.setstate() of version %s" %
142 (version, self.VERSION))
143
144## ---- 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
155 def __reduce__(self):
156 return self.__class__, (), self.getstate()
157
158## -------------------- integer methods -------------------
159
160 def randrange(self, start, stop=None, step=1, int=int, default=None,
161 maxwidth=1L<<BPF):
162 """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.
166 Do not supply the 'int', 'default', and 'maxwidth' arguments.
167 """
168
169 # This code is a bit messy to make it fast for the
170 # common case while still doing adequate error checking.
171 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:
176 if istart >= maxwidth:
177 return self._randbelow(istart)
178 return int(self.random() * istart)
179 raise ValueError, "empty range for randrange()"
180
181 # stop argument supplied.
182 istop = int(stop)
183 if istop != stop:
184 raise ValueError, "non-integer stop for randrange()"
185 width = istop - istart
186 if step == 1 and width > 0:
187 # Note that
188 # int(istart + self.random()*width)
189 # 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).
194 # istart + int(self.random()*width)
195 # 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).
199
200 if width >= maxwidth:
201 return int(istart + self._randbelow(width))
202 return int(istart + int(self.random()*width))
203 if step == 1:
204 raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width)
205
206 # Non-unit step argument supplied.
207 istep = int(step)
208 if istep != step:
209 raise ValueError, "non-integer step for randrange()"
210 if istep > 0:
211 n = (width + istep - 1) // istep
212 elif istep < 0:
213 n = (width + istep + 1) // istep
214 else:
215 raise ValueError, "zero step for randrange()"
216
217 if n <= 0:
218 raise ValueError, "empty range for randrange()"
219
220 if n >= maxwidth:
221 return istart + istep*self._randbelow(n)
222 return istart + istep*int(self.random() * n)
223
224 def randint(self, a, b):
225 """Return random integer in range [a, b], including both end points.
226 """
227
228 return self.randrange(a, b+1)
229
230 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
257## -------------------- sequence methods -------------------
258
259 def choice(self, seq):
260 """Choose a random element from a non-empty sequence."""
261 return seq[int(self.random() * len(seq))] # raises IndexError if seq is empty
262
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.
268 """
269
270 if random is None:
271 random = self.random
272 for i in reversed(xrange(1, len(x))):
273 # pick an element in x[:i+1] with which to exchange x[i]
274 j = int(random() * (i+1))
275 x[i], x[j] = x[j], x[i]
276
277 def sample(self, population, k):
278 """Chooses k unique random elements from a population sequence.
279
280 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).
285
286 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.
289
290 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)
293 """
294
295 # 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
304 # Sampling without replacement entails tracking either potential
305 # selections (the pool) in a list or previous selections in a set.
306
307 # When the number of selections is small compared to the
308 # population, then tracking selections is efficient, requiring
309 # only a small set and an occasional reselection. For
310 # a larger number of selections, the pool tracking method is
311 # preferred since the list takes less space than the
312 # set and it doesn't suffer from frequent reselections.
313
314 n = len(population)
315 if not 0 <= k <= n:
316 raise ValueError, "sample larger than population"
317 random = self.random
318 _int = int
319 result = [None] * k
320 setsize = 21 # size of a small set minus size of an empty list
321 if k > 5:
322 setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets
323 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.
326 pool = list(population)
327 for i in xrange(k): # invariant: non-selected at [0,n-i)
328 j = _int(random() * (n-i))
329 result[i] = pool[j]
330 pool[j] = pool[n-i-1] # move non-selected item into vacancy
331 else:
332 try:
333 selected = set()
334 selected_add = selected.add
335 for i in xrange(k):
336 j = _int(random() * n)
337 while j in selected:
338 j = _int(random() * n)
339 selected_add(j)
340 result[i] = population[j]
341 except (TypeError, KeyError): # handle (at least) sets
342 if isinstance(population, list):
343 raise
344 return self.sample(tuple(population), k)
345 return result
346
347## -------------------- real-valued distributions -------------------
348
349## -------------------- uniform distribution -------------------
350
351 def uniform(self, a, b):
352 "Get a random number in the range [a, b) or [a, b] depending on rounding."
353 return a + (b-a) * self.random()
354
355## -------------------- triangular --------------------
356
357 def triangular(self, low=0.0, high=1.0, mode=None):
358 """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()
367 c = 0.5 if mode is None else (mode - low) / (high - low)
368 if u > c:
369 u = 1.0 - u
370 c = 1.0 - c
371 low, high = high, low
372 return low + (high - low) * (u * c) ** 0.5
373
374## -------------------- normal distribution --------------------
375
376 def normalvariate(self, mu, sigma):
377 """Normal distribution.
378
379 mu is the mean, and sigma is the standard deviation.
380
381 """
382 # mu = mean, sigma = standard deviation
383
384 # 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.
388
389 random = self.random
390 while 1:
391 u1 = random()
392 u2 = 1.0 - random()
393 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
398
399## -------------------- lognormal distribution --------------------
400
401 def lognormvariate(self, mu, sigma):
402 """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.
407
408 """
409 return _exp(self.normalvariate(mu, sigma))
410
411## -------------------- exponential distribution --------------------
412
413 def expovariate(self, lambd):
414 """Exponential distribution.
415
416 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.
421
422 """
423 # lambd: rate lambd = 1/mean
424 # ('lambda' is a Python reserved word)
425
426 random = self.random
427 u = random()
428 while u <= 1e-7:
429 u = random()
430 return -_log(u)/lambd
431
432## -------------------- von Mises distribution --------------------
433
434 def vonmisesvariate(self, mu, kappa):
435 """Circular data distribution.
436
437 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.
441
442 """
443 # 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)
461
462 while 1:
463 u1 = random()
464
465 z = _cos(_pi * u1)
466 f = (1.0 + r * z)/(r + z)
467 c = kappa * (r - f)
468
469 u2 = random()
470
471 if u2 < c * (2.0 - c) or u2 <= c * _exp(1.0 - c):
472 break
473
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
482## -------------------- gamma distribution --------------------
483
484 def gammavariate(self, alpha, beta):
485 """Gamma distribution. Not the gamma function!
486
487 Conditions on the parameters are alpha > 0 and beta > 0.
488
489 """
490
491 # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
492
493 # 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'
497
498 random = self.random
499 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
505 ainv = _sqrt(2.0 * alpha - 1.0)
506 bbb = alpha - LOG4
507 ccc = alpha + ainv
508
509 while 1:
510 u1 = random()
511 if not 1e-7 < u1 < .9999999:
512 continue
513 u2 = 1.0 - random()
514 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):
519 return x * beta
520
521 elif alpha == 1.0:
522 # expovariate(1)
523 u = random()
524 while u <= 1e-7:
525 u = random()
526 return -_log(u) * beta
527
528 else: # alpha is between 0 and 1 (exclusive)
529
530 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
531
532 while 1:
533 u = random()
534 b = (_e + alpha)/_e
535 p = b*u
536 if p <= 1.0:
537 x = p ** (1.0/alpha)
538 else:
539 x = -_log((b-p)/alpha)
540 u1 = random()
541 if p > 1.0:
542 if u1 <= x ** (alpha - 1.0):
543 break
544 elif u1 <= _exp(-x):
545 break
546 return x * beta
547
548## -------------------- Gauss (faster alternative) --------------------
549
550 def gauss(self, mu, sigma):
551 """Gaussian distribution.
552
553 mu is the mean, and sigma is the standard deviation. This is
554 slightly faster than the normalvariate() function.
555
556 Not thread-safe without a lock around calls.
557
558 """
559
560 # 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)
570
571 # 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.)
577
578 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
586
587 return mu + z*sigma
588
589## -------------------- beta --------------------
590## 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.
602
603 def betavariate(self, alpha, beta):
604 """Beta distribution.
605
606 Conditions on the parameters are alpha > 0 and beta > 0.
607 Returned values range between 0 and 1.
608
609 """
610
611 # 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.))
618
619## -------------------- Pareto --------------------
620
621 def paretovariate(self, alpha):
622 """Pareto distribution. alpha is the shape parameter."""
623 # Jain, pg. 495
624
625 u = 1.0 - self.random()
626 return 1.0 / pow(u, 1.0/alpha)
627
628## -------------------- Weibull --------------------
629
630 def weibullvariate(self, alpha, beta):
631 """Weibull distribution.
632
633 alpha is the scale parameter and beta is the shape parameter.
634
635 """
636 # Jain, pg. 499; bug fix courtesy Bill Arms
637
638 u = 1.0 - self.random()
639 return alpha * pow(-_log(u), 1.0/beta)
640
641## -------------------- 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
650 None or no argument seeds from current time or from an operating
651 system specific randomness source if available.
652
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:
662 try:
663 a = long(_hexlify(_urandom(16)), 16)
664 except NotImplementedError:
665 import time
666 a = long(time.time() * 256) # use fractional seconds
667
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
791## --------------- Operating System Random Source ------------------
792
793class 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).
797
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)."""
803 return (long(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF
804
805 def getrandbits(self, k):
806 """getrandbits(k) -> x. Generates a long int with k random bits."""
807 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):
816 "Stub method. Not used for a system random number generator."
817 return None
818 seed = jumpahead = _stub
819
820 def _notimplemented(self, *args, **kwds):
821 "Method should not be called for a system random number generator."
822 raise NotImplementedError('System entropy source does not have state.')
823 getstate = setstate = _notimplemented
824
825## -------------------- test program --------------------
826
827def _test_generator(n, func, args):
828 import time
829 print n, 'times', func.__name__
830 total = 0.0
831 sqsum = 0.0
832 smallest = 1e10
833 largest = -1e10
834 t0 = time.time()
835 for i in range(n):
836 x = func(*args)
837 total += x
838 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,',
843 avg = total/n
844 stddev = _sqrt(sqsum/n - avg*avg)
845 print 'avg %g, stddev %g, min %g, max %g' % \
846 (avg, stddev, smallest, largest)
847
848
849def _test(N=2000):
850 _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))
865 _test_generator(N, triangular, (0.0, 1.0, 1.0/3.0))
866
867# Create one instance, seeded from current time, and export its methods
868# 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
873_inst = Random()
874seed = _inst.seed
875random = _inst.random
876uniform = _inst.uniform
877triangular = _inst.triangular
878randint = _inst.randint
879choice = _inst.choice
880randrange = _inst.randrange
881sample = _inst.sample
882shuffle = _inst.shuffle
883normalvariate = _inst.normalvariate
884lognormvariate = _inst.lognormvariate
885expovariate = _inst.expovariate
886vonmisesvariate = _inst.vonmisesvariate
887gammavariate = _inst.gammavariate
888gauss = _inst.gauss
889betavariate = _inst.betavariate
890paretovariate = _inst.paretovariate
891weibullvariate = _inst.weibullvariate
892getstate = _inst.getstate
893setstate = _inst.setstate
894jumpahead = _inst.jumpahead
895getrandbits = _inst.getrandbits
896
897if __name__ == '__main__':
898 _test()
Note: See TracBrowser for help on using the repository browser.