symbian-qemu-0.9.1-12/python-win32-2.6.1/lib/random.py
changeset 1 2fb8b9db1c86
equal deleted inserted replaced
0:ffa851df0825 1:2fb8b9db1c86
       
     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 
       
    30 General 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 
       
    42 from __future__ import division
       
    43 from warnings import warn as _warn
       
    44 from types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType
       
    45 from math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil
       
    46 from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
       
    47 from os import urandom as _urandom
       
    48 from 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 
       
    57 NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
       
    58 TWOPI = 2.0*_pi
       
    59 LOG4 = _log(4.0)
       
    60 SG_MAGICCONST = 1.0 + _log(4.5)
       
    61 BPF = 53        # Number of bits in a float
       
    62 RECIP_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 
       
    69 import _random
       
    70 
       
    71 class 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)."""
       
   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.  (The parameter would be
       
   417         called "lambda", but that is a reserved word in Python.)  Returned
       
   418         values range from 0 to positive infinity.
       
   419 
       
   420         """
       
   421         # lambd: rate lambd = 1/mean
       
   422         # ('lambda' is a Python reserved word)
       
   423 
       
   424         random = self.random
       
   425         u = random()
       
   426         while u <= 1e-7:
       
   427             u = random()
       
   428         return -_log(u)/lambd
       
   429 
       
   430 ## -------------------- von Mises distribution --------------------
       
   431 
       
   432     def vonmisesvariate(self, mu, kappa):
       
   433         """Circular data distribution.
       
   434 
       
   435         mu is the mean angle, expressed in radians between 0 and 2*pi, and
       
   436         kappa is the concentration parameter, which must be greater than or
       
   437         equal to zero.  If kappa is equal to zero, this distribution reduces
       
   438         to a uniform random angle over the range 0 to 2*pi.
       
   439 
       
   440         """
       
   441         # mu:    mean angle (in radians between 0 and 2*pi)
       
   442         # kappa: concentration parameter kappa (>= 0)
       
   443         # if kappa = 0 generate uniform random angle
       
   444 
       
   445         # Based upon an algorithm published in: Fisher, N.I.,
       
   446         # "Statistical Analysis of Circular Data", Cambridge
       
   447         # University Press, 1993.
       
   448 
       
   449         # Thanks to Magnus Kessler for a correction to the
       
   450         # implementation of step 4.
       
   451 
       
   452         random = self.random
       
   453         if kappa <= 1e-6:
       
   454             return TWOPI * random()
       
   455 
       
   456         a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa)
       
   457         b = (a - _sqrt(2.0 * a))/(2.0 * kappa)
       
   458         r = (1.0 + b * b)/(2.0 * b)
       
   459 
       
   460         while 1:
       
   461             u1 = random()
       
   462 
       
   463             z = _cos(_pi * u1)
       
   464             f = (1.0 + r * z)/(r + z)
       
   465             c = kappa * (r - f)
       
   466 
       
   467             u2 = random()
       
   468 
       
   469             if u2 < c * (2.0 - c) or u2 <= c * _exp(1.0 - c):
       
   470                 break
       
   471 
       
   472         u3 = random()
       
   473         if u3 > 0.5:
       
   474             theta = (mu % TWOPI) + _acos(f)
       
   475         else:
       
   476             theta = (mu % TWOPI) - _acos(f)
       
   477 
       
   478         return theta
       
   479 
       
   480 ## -------------------- gamma distribution --------------------
       
   481 
       
   482     def gammavariate(self, alpha, beta):
       
   483         """Gamma distribution.  Not the gamma function!
       
   484 
       
   485         Conditions on the parameters are alpha > 0 and beta > 0.
       
   486 
       
   487         """
       
   488 
       
   489         # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
       
   490 
       
   491         # Warning: a few older sources define the gamma distribution in terms
       
   492         # of alpha > -1.0
       
   493         if alpha <= 0.0 or beta <= 0.0:
       
   494             raise ValueError, 'gammavariate: alpha and beta must be > 0.0'
       
   495 
       
   496         random = self.random
       
   497         if alpha > 1.0:
       
   498 
       
   499             # Uses R.C.H. Cheng, "The generation of Gamma
       
   500             # variables with non-integral shape parameters",
       
   501             # Applied Statistics, (1977), 26, No. 1, p71-74
       
   502 
       
   503             ainv = _sqrt(2.0 * alpha - 1.0)
       
   504             bbb = alpha - LOG4
       
   505             ccc = alpha + ainv
       
   506 
       
   507             while 1:
       
   508                 u1 = random()
       
   509                 if not 1e-7 < u1 < .9999999:
       
   510                     continue
       
   511                 u2 = 1.0 - random()
       
   512                 v = _log(u1/(1.0-u1))/ainv
       
   513                 x = alpha*_exp(v)
       
   514                 z = u1*u1*u2
       
   515                 r = bbb+ccc*v-x
       
   516                 if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z):
       
   517                     return x * beta
       
   518 
       
   519         elif alpha == 1.0:
       
   520             # expovariate(1)
       
   521             u = random()
       
   522             while u <= 1e-7:
       
   523                 u = random()
       
   524             return -_log(u) * beta
       
   525 
       
   526         else:   # alpha is between 0 and 1 (exclusive)
       
   527 
       
   528             # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
       
   529 
       
   530             while 1:
       
   531                 u = random()
       
   532                 b = (_e + alpha)/_e
       
   533                 p = b*u
       
   534                 if p <= 1.0:
       
   535                     x = p ** (1.0/alpha)
       
   536                 else:
       
   537                     x = -_log((b-p)/alpha)
       
   538                 u1 = random()
       
   539                 if p > 1.0:
       
   540                     if u1 <= x ** (alpha - 1.0):
       
   541                         break
       
   542                 elif u1 <= _exp(-x):
       
   543                     break
       
   544             return x * beta
       
   545 
       
   546 ## -------------------- Gauss (faster alternative) --------------------
       
   547 
       
   548     def gauss(self, mu, sigma):
       
   549         """Gaussian distribution.
       
   550 
       
   551         mu is the mean, and sigma is the standard deviation.  This is
       
   552         slightly faster than the normalvariate() function.
       
   553 
       
   554         Not thread-safe without a lock around calls.
       
   555 
       
   556         """
       
   557 
       
   558         # When x and y are two variables from [0, 1), uniformly
       
   559         # distributed, then
       
   560         #
       
   561         #    cos(2*pi*x)*sqrt(-2*log(1-y))
       
   562         #    sin(2*pi*x)*sqrt(-2*log(1-y))
       
   563         #
       
   564         # are two *independent* variables with normal distribution
       
   565         # (mu = 0, sigma = 1).
       
   566         # (Lambert Meertens)
       
   567         # (corrected version; bug discovered by Mike Miller, fixed by LM)
       
   568 
       
   569         # Multithreading note: When two threads call this function
       
   570         # simultaneously, it is possible that they will receive the
       
   571         # same return value.  The window is very small though.  To
       
   572         # avoid this, you have to use a lock around all calls.  (I
       
   573         # didn't want to slow this down in the serial case by using a
       
   574         # lock here.)
       
   575 
       
   576         random = self.random
       
   577         z = self.gauss_next
       
   578         self.gauss_next = None
       
   579         if z is None:
       
   580             x2pi = random() * TWOPI
       
   581             g2rad = _sqrt(-2.0 * _log(1.0 - random()))
       
   582             z = _cos(x2pi) * g2rad
       
   583             self.gauss_next = _sin(x2pi) * g2rad
       
   584 
       
   585         return mu + z*sigma
       
   586 
       
   587 ## -------------------- beta --------------------
       
   588 ## See
       
   589 ## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470
       
   590 ## for Ivan Frohne's insightful analysis of why the original implementation:
       
   591 ##
       
   592 ##    def betavariate(self, alpha, beta):
       
   593 ##        # Discrete Event Simulation in C, pp 87-88.
       
   594 ##
       
   595 ##        y = self.expovariate(alpha)
       
   596 ##        z = self.expovariate(1.0/beta)
       
   597 ##        return z/(y+z)
       
   598 ##
       
   599 ## was dead wrong, and how it probably got that way.
       
   600 
       
   601     def betavariate(self, alpha, beta):
       
   602         """Beta distribution.
       
   603 
       
   604         Conditions on the parameters are alpha > 0 and beta > 0.
       
   605         Returned values range between 0 and 1.
       
   606 
       
   607         """
       
   608 
       
   609         # This version due to Janne Sinkkonen, and matches all the std
       
   610         # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
       
   611         y = self.gammavariate(alpha, 1.)
       
   612         if y == 0:
       
   613             return 0.0
       
   614         else:
       
   615             return y / (y + self.gammavariate(beta, 1.))
       
   616 
       
   617 ## -------------------- Pareto --------------------
       
   618 
       
   619     def paretovariate(self, alpha):
       
   620         """Pareto distribution.  alpha is the shape parameter."""
       
   621         # Jain, pg. 495
       
   622 
       
   623         u = 1.0 - self.random()
       
   624         return 1.0 / pow(u, 1.0/alpha)
       
   625 
       
   626 ## -------------------- Weibull --------------------
       
   627 
       
   628     def weibullvariate(self, alpha, beta):
       
   629         """Weibull distribution.
       
   630 
       
   631         alpha is the scale parameter and beta is the shape parameter.
       
   632 
       
   633         """
       
   634         # Jain, pg. 499; bug fix courtesy Bill Arms
       
   635 
       
   636         u = 1.0 - self.random()
       
   637         return alpha * pow(-_log(u), 1.0/beta)
       
   638 
       
   639 ## -------------------- Wichmann-Hill -------------------
       
   640 
       
   641 class WichmannHill(Random):
       
   642 
       
   643     VERSION = 1     # used by getstate/setstate
       
   644 
       
   645     def seed(self, a=None):
       
   646         """Initialize internal state from hashable object.
       
   647 
       
   648         None or no argument seeds from current time or from an operating
       
   649         system specific randomness source if available.
       
   650 
       
   651         If a is not None or an int or long, hash(a) is used instead.
       
   652 
       
   653         If a is an int or long, a is used directly.  Distinct values between
       
   654         0 and 27814431486575L inclusive are guaranteed to yield distinct
       
   655         internal states (this guarantee is specific to the default
       
   656         Wichmann-Hill generator).
       
   657         """
       
   658 
       
   659         if a is None:
       
   660             try:
       
   661                 a = long(_hexlify(_urandom(16)), 16)
       
   662             except NotImplementedError:
       
   663                 import time
       
   664                 a = long(time.time() * 256) # use fractional seconds
       
   665 
       
   666         if not isinstance(a, (int, long)):
       
   667             a = hash(a)
       
   668 
       
   669         a, x = divmod(a, 30268)
       
   670         a, y = divmod(a, 30306)
       
   671         a, z = divmod(a, 30322)
       
   672         self._seed = int(x)+1, int(y)+1, int(z)+1
       
   673 
       
   674         self.gauss_next = None
       
   675 
       
   676     def random(self):
       
   677         """Get the next random number in the range [0.0, 1.0)."""
       
   678 
       
   679         # Wichman-Hill random number generator.
       
   680         #
       
   681         # Wichmann, B. A. & Hill, I. D. (1982)
       
   682         # Algorithm AS 183:
       
   683         # An efficient and portable pseudo-random number generator
       
   684         # Applied Statistics 31 (1982) 188-190
       
   685         #
       
   686         # see also:
       
   687         #        Correction to Algorithm AS 183
       
   688         #        Applied Statistics 33 (1984) 123
       
   689         #
       
   690         #        McLeod, A. I. (1985)
       
   691         #        A remark on Algorithm AS 183
       
   692         #        Applied Statistics 34 (1985),198-200
       
   693 
       
   694         # This part is thread-unsafe:
       
   695         # BEGIN CRITICAL SECTION
       
   696         x, y, z = self._seed
       
   697         x = (171 * x) % 30269
       
   698         y = (172 * y) % 30307
       
   699         z = (170 * z) % 30323
       
   700         self._seed = x, y, z
       
   701         # END CRITICAL SECTION
       
   702 
       
   703         # Note:  on a platform using IEEE-754 double arithmetic, this can
       
   704         # never return 0.0 (asserted by Tim; proof too long for a comment).
       
   705         return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0
       
   706 
       
   707     def getstate(self):
       
   708         """Return internal state; can be passed to setstate() later."""
       
   709         return self.VERSION, self._seed, self.gauss_next
       
   710 
       
   711     def setstate(self, state):
       
   712         """Restore internal state from object returned by getstate()."""
       
   713         version = state[0]
       
   714         if version == 1:
       
   715             version, self._seed, self.gauss_next = state
       
   716         else:
       
   717             raise ValueError("state with version %s passed to "
       
   718                              "Random.setstate() of version %s" %
       
   719                              (version, self.VERSION))
       
   720 
       
   721     def jumpahead(self, n):
       
   722         """Act as if n calls to random() were made, but quickly.
       
   723 
       
   724         n is an int, greater than or equal to 0.
       
   725 
       
   726         Example use:  If you have 2 threads and know that each will
       
   727         consume no more than a million random numbers, create two Random
       
   728         objects r1 and r2, then do
       
   729             r2.setstate(r1.getstate())
       
   730             r2.jumpahead(1000000)
       
   731         Then r1 and r2 will use guaranteed-disjoint segments of the full
       
   732         period.
       
   733         """
       
   734 
       
   735         if not n >= 0:
       
   736             raise ValueError("n must be >= 0")
       
   737         x, y, z = self._seed
       
   738         x = int(x * pow(171, n, 30269)) % 30269
       
   739         y = int(y * pow(172, n, 30307)) % 30307
       
   740         z = int(z * pow(170, n, 30323)) % 30323
       
   741         self._seed = x, y, z
       
   742 
       
   743     def __whseed(self, x=0, y=0, z=0):
       
   744         """Set the Wichmann-Hill seed from (x, y, z).
       
   745 
       
   746         These must be integers in the range [0, 256).
       
   747         """
       
   748 
       
   749         if not type(x) == type(y) == type(z) == int:
       
   750             raise TypeError('seeds must be integers')
       
   751         if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256):
       
   752             raise ValueError('seeds must be in range(0, 256)')
       
   753         if 0 == x == y == z:
       
   754             # Initialize from current time
       
   755             import time
       
   756             t = long(time.time() * 256)
       
   757             t = int((t&0xffffff) ^ (t>>24))
       
   758             t, x = divmod(t, 256)
       
   759             t, y = divmod(t, 256)
       
   760             t, z = divmod(t, 256)
       
   761         # Zero is a poor seed, so substitute 1
       
   762         self._seed = (x or 1, y or 1, z or 1)
       
   763 
       
   764         self.gauss_next = None
       
   765 
       
   766     def whseed(self, a=None):
       
   767         """Seed from hashable object's hash code.
       
   768 
       
   769         None or no argument seeds from current time.  It is not guaranteed
       
   770         that objects with distinct hash codes lead to distinct internal
       
   771         states.
       
   772 
       
   773         This is obsolete, provided for compatibility with the seed routine
       
   774         used prior to Python 2.1.  Use the .seed() method instead.
       
   775         """
       
   776 
       
   777         if a is None:
       
   778             self.__whseed()
       
   779             return
       
   780         a = hash(a)
       
   781         a, x = divmod(a, 256)
       
   782         a, y = divmod(a, 256)
       
   783         a, z = divmod(a, 256)
       
   784         x = (x + a) % 256 or 1
       
   785         y = (y + a) % 256 or 1
       
   786         z = (z + a) % 256 or 1
       
   787         self.__whseed(x, y, z)
       
   788 
       
   789 ## --------------- Operating System Random Source  ------------------
       
   790 
       
   791 class SystemRandom(Random):
       
   792     """Alternate random number generator using sources provided
       
   793     by the operating system (such as /dev/urandom on Unix or
       
   794     CryptGenRandom on Windows).
       
   795 
       
   796      Not available on all systems (see os.urandom() for details).
       
   797     """
       
   798 
       
   799     def random(self):
       
   800         """Get the next random number in the range [0.0, 1.0)."""
       
   801         return (long(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF
       
   802 
       
   803     def getrandbits(self, k):
       
   804         """getrandbits(k) -> x.  Generates a long int with k random bits."""
       
   805         if k <= 0:
       
   806             raise ValueError('number of bits must be greater than zero')
       
   807         if k != int(k):
       
   808             raise TypeError('number of bits should be an integer')
       
   809         bytes = (k + 7) // 8                    # bits / 8 and rounded up
       
   810         x = long(_hexlify(_urandom(bytes)), 16)
       
   811         return x >> (bytes * 8 - k)             # trim excess bits
       
   812 
       
   813     def _stub(self, *args, **kwds):
       
   814         "Stub method.  Not used for a system random number generator."
       
   815         return None
       
   816     seed = jumpahead = _stub
       
   817 
       
   818     def _notimplemented(self, *args, **kwds):
       
   819         "Method should not be called for a system random number generator."
       
   820         raise NotImplementedError('System entropy source does not have state.')
       
   821     getstate = setstate = _notimplemented
       
   822 
       
   823 ## -------------------- test program --------------------
       
   824 
       
   825 def _test_generator(n, func, args):
       
   826     import time
       
   827     print n, 'times', func.__name__
       
   828     total = 0.0
       
   829     sqsum = 0.0
       
   830     smallest = 1e10
       
   831     largest = -1e10
       
   832     t0 = time.time()
       
   833     for i in range(n):
       
   834         x = func(*args)
       
   835         total += x
       
   836         sqsum = sqsum + x*x
       
   837         smallest = min(x, smallest)
       
   838         largest = max(x, largest)
       
   839     t1 = time.time()
       
   840     print round(t1-t0, 3), 'sec,',
       
   841     avg = total/n
       
   842     stddev = _sqrt(sqsum/n - avg*avg)
       
   843     print 'avg %g, stddev %g, min %g, max %g' % \
       
   844               (avg, stddev, smallest, largest)
       
   845 
       
   846 
       
   847 def _test(N=2000):
       
   848     _test_generator(N, random, ())
       
   849     _test_generator(N, normalvariate, (0.0, 1.0))
       
   850     _test_generator(N, lognormvariate, (0.0, 1.0))
       
   851     _test_generator(N, vonmisesvariate, (0.0, 1.0))
       
   852     _test_generator(N, gammavariate, (0.01, 1.0))
       
   853     _test_generator(N, gammavariate, (0.1, 1.0))
       
   854     _test_generator(N, gammavariate, (0.1, 2.0))
       
   855     _test_generator(N, gammavariate, (0.5, 1.0))
       
   856     _test_generator(N, gammavariate, (0.9, 1.0))
       
   857     _test_generator(N, gammavariate, (1.0, 1.0))
       
   858     _test_generator(N, gammavariate, (2.0, 1.0))
       
   859     _test_generator(N, gammavariate, (20.0, 1.0))
       
   860     _test_generator(N, gammavariate, (200.0, 1.0))
       
   861     _test_generator(N, gauss, (0.0, 1.0))
       
   862     _test_generator(N, betavariate, (3.0, 3.0))
       
   863     _test_generator(N, triangular, (0.0, 1.0, 1.0/3.0))
       
   864 
       
   865 # Create one instance, seeded from current time, and export its methods
       
   866 # as module-level functions.  The functions share state across all uses
       
   867 #(both in the user's code and in the Python libraries), but that's fine
       
   868 # for most programs and is easier for the casual user than making them
       
   869 # instantiate their own Random() instance.
       
   870 
       
   871 _inst = Random()
       
   872 seed = _inst.seed
       
   873 random = _inst.random
       
   874 uniform = _inst.uniform
       
   875 triangular = _inst.triangular
       
   876 randint = _inst.randint
       
   877 choice = _inst.choice
       
   878 randrange = _inst.randrange
       
   879 sample = _inst.sample
       
   880 shuffle = _inst.shuffle
       
   881 normalvariate = _inst.normalvariate
       
   882 lognormvariate = _inst.lognormvariate
       
   883 expovariate = _inst.expovariate
       
   884 vonmisesvariate = _inst.vonmisesvariate
       
   885 gammavariate = _inst.gammavariate
       
   886 gauss = _inst.gauss
       
   887 betavariate = _inst.betavariate
       
   888 paretovariate = _inst.paretovariate
       
   889 weibullvariate = _inst.weibullvariate
       
   890 getstate = _inst.getstate
       
   891 setstate = _inst.setstate
       
   892 jumpahead = _inst.jumpahead
       
   893 getrandbits = _inst.getrandbits
       
   894 
       
   895 if __name__ == '__main__':
       
   896     _test()