Reworked random.py so that it no longer depends on, and offers all the

functionality of, whrandom.py.  Also closes all the "XXX" todos in
random.py.  New frequently-requested functions/methods getstate() and
setstate().  All exported functions are now bound methods of a hidden
instance.  Killed all unintended exports.  Updated the docs.
FRED:  The more I fiddle the docs, the less I understand the exact
intended use of the \var, \code, \method tags.  Please review critically.
GUIDO:  See email.  I updated NEWS as if whrandom were deprecated; I
think it should be.

3 files changed, 481 insertions, 288 deletions

diff --git a/Doc/lib/librandom.tex b/Doc/lib/librandom.tex
index 82a55e6..4e4d615 100644
--- a/Doc/lib/librandom.tex
+++ b/Doc/lib/librandom.tex
@@ -32,6 +32,18 @@ functions from multiple threads, you should explicitly serialize the calls.
 Else, because no critical sections are implemented internally, calls
 from different threads may see the same return values.
 
+The functions supplied by this module are actually bound methods of a
+hidden instance of the \var{random.Random} class.  You can instantiate
+your own instances of \var{Random} to get generators that don't share state.
+This may be especially useful for multi-threaded programs, although there's
+no simple way to seed the distinct generators to ensure that the generated
+sequences won't overlap.  Class \var{Random} can also be subclassed if you
+want to use a different basic generator of your own devising:  in that
+case, override the \method{random()}, \method{seed()}, \method{getstate()}
+and \method{setstate()} methods.
+
+
+Bookkeeping functions:
 
 \begin{funcdesc}{seed}{\optional{x}}
   Initialize the basic random number generator.
@@ -45,15 +57,19 @@ from different threads may see the same return values.
   the module is first imported.
 \end{funcdesc}
 
-\begin{funcdesc}{choice}{seq}
-  Return a random element from the non-empty sequence \var{seq}.
-\end{funcdesc}
+\begin{funcdesc}{getstate}{}
+  Return an object capturing the current internal state of the generator.
+  This object can be passed to \code{setstate()} to restore the state.
+ \end{funcdesc}
+
+\begin{funcdesc}{setstate}{state}
+  \var{state} should have been obtained from a previous call to
+  \code{getstate()}, and \code{setstate()} restores the internal state
+  of the generate to what it was at the time \code{setstate()} was called.
+ \end{funcdesc}
 
-\begin{funcdesc}{randint}{a, b}
-  \deprecated{2.0}{Use \function{randrange()} instead.}
-  Return a random integer \var{N} such that
-  \code{\var{a} <= \var{N} <= \var{b}}.
-\end{funcdesc}
+
+Functions for integers:
 
 \begin{funcdesc}{randrange}{\optional{start,} stop\optional{, step}}
   Return a randomly selected element from \code{range(\var{start},
@@ -63,6 +79,37 @@ from different threads may see the same return values.
   \versionadded{1.5.2}
 \end{funcdesc}
 
+\begin{funcdesc}{randint}{a, b}
+  \deprecated{2.0}{Use \function{randrange()} instead.}
+  Return a random integer \var{N} such that
+  \code{\var{a} <= \var{N} <= \var{b}}.
+\end{funcdesc}
+
+
+Functions for sequences:
+
+\begin{funcdesc}{choice}{seq}
+  Return a random element from the non-empty sequence \var{seq}.
+\end{funcdesc}
+
+\begin{funcdesc}{shuffle}{x\optional{, random}}
+  Shuffle the sequence \var{x} in place.
+  The optional argument \var{random} is a 0-argument function
+  returning a random float in [0.0, 1.0); by default, this is the
+  function \function{random()}.
+
+  Note that for even rather small \code{len(\var{x})}, the total
+  number of permutations of \var{x} is larger than the period of most
+  random number generators; this implies that most permutations of a
+  long sequence can never be generated.
+\end{funcdesc}
+
+
+The following functions generate specific real-valued distributions.
+Function parameters are named after the corresponding variables in the
+distribution's equation, as used in common mathematical practice; most of
+these equations can be found in any statistics text.
+
 \begin{funcdesc}{random}{}
   Return the next random floating point number in the range [0.0, 1.0).
 \end{funcdesc}
@@ -72,14 +119,6 @@ from different threads may see the same return values.
   \code{\var{a} <= \var{N} < \var{b}}.
 \end{funcdesc}
 
-
-The following functions are defined to support specific distributions,
-and all return real values.  Function parameters are named after the
-corresponding variables in the distribution's equation, as used in
-common mathematical practice; most of these equations can be found in
-any statistics text.
-
-
 \begin{funcdesc}{betavariate}{alpha, beta}
   Beta distribution.  Conditions on the parameters are
   \code{\var{alpha} > -1} and \code{\var{beta} > -1}.
@@ -143,21 +182,6 @@ any statistics text.
 \end{funcdesc}
 
 
-This function does not represent a specific distribution, but
-implements a standard useful algorithm:
-
-\begin{funcdesc}{shuffle}{x\optional{, random}}
-  Shuffle the sequence \var{x} in place.
-  The optional argument \var{random} is a 0-argument function
-  returning a random float in [0.0, 1.0); by default, this is the
-  function \function{random()}.
-
-  Note that for even rather small \code{len(\var{x})}, the total
-  number of permutations of \var{x} is larger than the period of most
-  random number generators; this implies that most permutations of a
-  long sequence can never be generated.
-\end{funcdesc}
-
 \begin{seealso}
   \seetext{Wichmann, B. A. \& Hill, I. D., ``Algorithm AS 183:
            An efficient and portable pseudo-random number generator'',
diff --git a/Lib/random.py b/Lib/random.py
index d10ce78..a818f73 100644
--- a/Lib/random.py
+++ b/Lib/random.py
@@ -1,7 +1,17 @@
 """Random variable generators.
 
+    integers
+    --------
+           uniform within range
+
+    sequences
+    ---------
+           pick random element
+           generate random permutation
+
     distributions on the real line:
     ------------------------------
+           uniform
            normal (Gaussian)
            lognormal
            negative exponential
@@ -17,328 +27,429 @@ Translated from anonymously contributed C/C++ source.
 
 Multi-threading note: the random number generator used here is not
 thread-safe; it is possible that two calls return the same random
-value.  See whrandom.py for more info.
+value.
 """
+# XXX The docstring sucks.
 
-import whrandom
-from whrandom import random, uniform, randint, choice, randrange # For export!
-from math import log, exp, pi, e, sqrt, acos, cos, sin
+from math import log as _log, exp as _exp, pi as _pi, e as _e
+from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
+
+def _verify(name, expected):
+    computed = eval(name)
+    if abs(computed - expected) > 1e-7:
+        raise ValueError(
+            "computed value for %s deviates too much "
+            "(computed %g, expected %g)" % (name, computed, expected))
 
-# Interfaces to replace remaining needs for importing whrandom
-# XXX TO DO: make the distribution functions below into methods.
+NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
+_verify('NV_MAGICCONST', 1.71552776992141)
 
-def makeseed(a=None):
-    """Turn a hashable value into three seed values for whrandom.seed().
+TWOPI = 2.0*_pi
+_verify('TWOPI', 6.28318530718)
 
-    None or no argument returns (0, 0, 0), to seed from current time.
+LOG4 = _log(4.0)
+_verify('LOG4', 1.38629436111989)
 
-    """
-    if a is None:
-        return (0, 0, 0)
-    a = hash(a)
-    a, x = divmod(a, 256)
-    a, y = divmod(a, 256)
-    a, z = divmod(a, 256)
-    x = (x + a) % 256 or 1
-    y = (y + a) % 256 or 1
-    z = (z + a) % 256 or 1
-    return (x, y, z)
+SG_MAGICCONST = 1.0 + _log(4.5)
+_verify('SG_MAGICCONST', 2.50407739677627)
 
-def seed(a=None):
-    """Seed the default generator from any hashable value.
+del _verify
 
-    None or no argument seeds from current time.
+# Translated by Guido van Rossum from C source provided by
+# Adrian Baddeley.
 
-    """
-    x, y, z = makeseed(a)
-    whrandom.seed(x, y, z)
+class Random:
 
-class generator(whrandom.whrandom):
-    """Random generator class."""
+    VERSION = 1     # used by getstate/setstate
 
-    def __init__(self, a=None):
-        """Constructor.  Seed from current time or hashable value."""
-        self.seed(a)
+    def __init__(self, x=None):
+        """Initialize an instance.
 
-    def seed(self, a=None):
-        """Seed the generator from current time or hashable value."""
-        x, y, z = makeseed(a)
-        whrandom.whrandom.seed(self, x, y, z)
+        Optional argument x controls seeding, as for Random.seed().
+        """
 
-def new_generator(a=None):
-    """Return a new random generator instance."""
-    return generator(a)
+        self.seed(x)
+        self.gauss_next = None
 
-# Housekeeping function to verify that magic constants have been
-# computed correctly
+    # Specific to Wichmann-Hill generator.  Subclasses wishing to use a
+    # different core generator should override seed(), random(),  getstate()
+    # and setstate().
 
-def verify(name, expected):
-    computed = eval(name)
-    if abs(computed - expected) > 1e-7:
-        raise ValueError, \
-'computed value for %s deviates too much (computed %g, expected %g)' % \
-(name, computed, expected)
+    def __whseed(self, x=0, y=0, z=0):
+        """Set the Wichmann-Hill seed from (x, y, z).
+
+        These must be integers in the range [0, 256).
+        """
+
+        if not type(x) == type(y) == type(z) == type(0):
+            raise TypeError('seeds must be integers')
+        if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256):
+            raise ValueError('seeds must be in range(0, 256)')
+        if 0 == x == y == z:
+            # Initialize from current time
+            import time
+            t = long(time.time()) * 256
+            t = int((t&0xffffff) ^ (t>>24))
+            t, x = divmod(t, 256)
+            t, y = divmod(t, 256)
+            t, z = divmod(t, 256)
+        # Zero is a poor seed, so substitute 1
+        self._seed = (x or 1, y or 1, z or 1)
+
+    def seed(self, a=None):
+        """Seed from hashable value
+
+        None or no argument seeds from current time.
+        """
+
+        if a is None:
+            self.__whseed()
+        a = hash(a)
+        a, x = divmod(a, 256)
+        a, y = divmod(a, 256)
+        a, z = divmod(a, 256)
+        x = (x + a) % 256 or 1
+        y = (y + a) % 256 or 1
+        z = (z + a) % 256 or 1
+        self.__whseed(x, y, z)
+
+    def getstate(self):
+        """Return internal state; can be passed to setstate() later."""
+
+        return self.VERSION, self._seed, self.gauss_next
+
+    def __getstate__(self): # for pickle
+        self.getstate()
+
+    def setstate(self, state):
+        """Restore internal state from object returned by getstate()."""
+        version = state[0]
+        if version == 1:
+            version, self._seed, self.gauss_next = state
+        else:
+            raise ValueError("state with version %s passed to "
+                             "Random.setstate() of version %s" %
+                             (version, self.VERSION))
+
+    def __setstate__(self, state):  # for pickle
+        self.setstate(state)
+
+    def random(self):
+        """Get the next random number in the range [0.0, 1.0)."""
+
+        # Wichman-Hill random number generator.
+        #
+        # Wichmann, B. A. & Hill, I. D. (1982)
+        # Algorithm AS 183:
+        # An efficient and portable pseudo-random number generator
+        # Applied Statistics 31 (1982) 188-190
+        #
+        # see also:
+        #        Correction to Algorithm AS 183
+        #        Applied Statistics 33 (1984) 123
+        #
+        #        McLeod, A. I. (1985)
+        #        A remark on Algorithm AS 183
+        #        Applied Statistics 34 (1985),198-200
+
+        # This part is thread-unsafe:
+        # BEGIN CRITICAL SECTION
+        x, y, z = self._seed
+        x = (171 * x) % 30269
+        y = (172 * y) % 30307
+        z = (170 * z) % 30323
+        self._seed = x, y, z
+        # END CRITICAL SECTION
+
+        # Note:  on a platform using IEEE-754 double arithmetic, this can
+        # never return 0.0 (asserted by Tim; proof too long for a comment).
+        return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0
+
+    def randrange(self, start, stop=None, step=1, int=int, default=None):
+        """Choose a random item from range(start, stop[, step]).
+
+        This fixes the problem with randint() which includes the
+        endpoint; in Python this is usually not what you want.
+        Do not supply the 'int' and 'default' arguments.
+        """
+
+        # This code is a bit messy to make it fast for the
+        # common case while still doing adequate error checking
+        istart = int(start)
+        if istart != start:
+            raise ValueError, "non-integer arg 1 for randrange()"
+        if stop is default:
+            if istart > 0:
+                return int(self.random() * istart)
+            raise ValueError, "empty range for randrange()"
+        istop = int(stop)
+        if istop != stop:
+            raise ValueError, "non-integer stop for randrange()"
+        if step == 1:
+            if istart < istop:
+                return istart + int(self.random() *
+                                   (istop - istart))
+            raise ValueError, "empty range for randrange()"
+        istep = int(step)
+        if istep != step:
+            raise ValueError, "non-integer step for randrange()"
+        if istep > 0:
+            n = (istop - istart + istep - 1) / istep
+        elif istep < 0:
+            n = (istop - istart + istep + 1) / istep
+        else:
+            raise ValueError, "zero step for randrange()"
+
+        if n <= 0:
+            raise ValueError, "empty range for randrange()"
+        return istart + istep*int(self.random() * n)
+
+    def randint(self, a, b):
+        """Get a random integer in the range [a, b] including
+        both end points.
+
+        (Deprecated; use randrange below.)
+        """
+
+        return self.randrange(a, b+1)
+
+    def choice(self, seq):
+        """Choose a random element from a non-empty sequence."""
+        return seq[int(self.random() * len(seq))]
+
+    def shuffle(self, x, random=None, int=int):
+        """x, random=random.random -> shuffle list x in place; return None.
+
+        Optional arg random is a 0-argument function returning a random
+        float in [0.0, 1.0); by default, the standard random.random.
+
+        Note that for even rather small len(x), the total number of
+        permutations of x is larger than the period of most random number
+        generators; this implies that "most" permutations of a long
+        sequence can never be generated.
+        """
+
+        if random is None:
+            random = self.random
+        for i in xrange(len(x)-1, 0, -1):
+        # pick an element in x[:i+1] with which to exchange x[i]
+            j = int(random() * (i+1))
+            x[i], x[j] = x[j], x[i]
+
+# -------------------- uniform distribution -------------------
+
+    def uniform(self, a, b):
+        """Get a random number in the range [a, b)."""
+        return a + (b-a) * self.random()
 
 # -------------------- normal distribution --------------------
 
-NV_MAGICCONST = 4*exp(-0.5)/sqrt(2.0)
-verify('NV_MAGICCONST', 1.71552776992141)
-def normalvariate(mu, sigma):
-    # mu = mean, sigma = standard deviation
-
-    # Uses Kinderman and Monahan method. Reference: Kinderman,
-    # A.J. and Monahan, J.F., "Computer generation of random
-    # variables using the ratio of uniform deviates", ACM Trans
-    # Math Software, 3, (1977), pp257-260.
-
-    while 1:
-        u1 = random()
-        u2 = random()
-        z = NV_MAGICCONST*(u1-0.5)/u2
-        zz = z*z/4.0
-        if zz <= -log(u2):
-            break
-    return mu+z*sigma
+    def normalvariate(self, mu, sigma):
+        # mu = mean, sigma = standard deviation
+
+        # Uses Kinderman and Monahan method. Reference: Kinderman,
+        # A.J. and Monahan, J.F., "Computer generation of random
+        # variables using the ratio of uniform deviates", ACM Trans
+        # Math Software, 3, (1977), pp257-260.
+
+        random = self.random
+        while 1:
+            u1 = random()
+            u2 = random()
+            z = NV_MAGICCONST*(u1-0.5)/u2
+            zz = z*z/4.0
+            if zz <= -_log(u2):
+                break
+        return mu + z*sigma
 
 # -------------------- lognormal distribution --------------------
 
-def lognormvariate(mu, sigma):
-    return exp(normalvariate(mu, sigma))
+    def lognormvariate(self, mu, sigma):
+        return _exp(self.normalvariate(mu, sigma))
 
 # -------------------- circular uniform --------------------
 
-def cunifvariate(mean, arc):
-    # mean: mean angle (in radians between 0 and pi)
-    # arc:  range of distribution (in radians between 0 and pi)
+    def cunifvariate(self, mean, arc):
+        # mean: mean angle (in radians between 0 and pi)
+        # arc:  range of distribution (in radians between 0 and pi)
 
-    return (mean + arc * (random() - 0.5)) % pi
+        return (mean + arc * (self.random() - 0.5)) % _pi
 
 # -------------------- exponential distribution --------------------
 
-def expovariate(lambd):
-    # lambd: rate lambd = 1/mean
-    # ('lambda' is a Python reserved word)
+    def expovariate(self, lambd):
+        # lambd: rate lambd = 1/mean
+        # ('lambda' is a Python reserved word)
 
-    u = random()
-    while u <= 1e-7:
+        random = self.random
         u = random()
-    return -log(u)/lambd
+        while u <= 1e-7:
+            u = random()
+        return -_log(u)/lambd
 
 # -------------------- von Mises distribution --------------------
 
-TWOPI = 2.0*pi
-verify('TWOPI', 6.28318530718)
-
-def vonmisesvariate(mu, kappa):
-    # mu:    mean angle (in radians between 0 and 2*pi)
-    # kappa: concentration parameter kappa (>= 0)
-    # if kappa = 0 generate uniform random angle
+    def vonmisesvariate(self, mu, kappa):
+        # mu:    mean angle (in radians between 0 and 2*pi)
+        # kappa: concentration parameter kappa (>= 0)
+        # if kappa = 0 generate uniform random angle
 
-    # Based upon an algorithm published in: Fisher, N.I.,
-    # "Statistical Analysis of Circular Data", Cambridge
-    # University Press, 1993.
+        # Based upon an algorithm published in: Fisher, N.I.,
+        # "Statistical Analysis of Circular Data", Cambridge
+        # University Press, 1993.
 
-    # Thanks to Magnus Kessler for a correction to the
-    # implementation of step 4.
+        # Thanks to Magnus Kessler for a correction to the
+        # implementation of step 4.
 
-    if kappa <= 1e-6:
-        return TWOPI * random()
+        random = self.random
+        if kappa <= 1e-6:
+            return TWOPI * random()
 
-    a = 1.0 + sqrt(1.0 + 4.0 * kappa * kappa)
-    b = (a - sqrt(2.0 * a))/(2.0 * kappa)
-    r = (1.0 + b * b)/(2.0 * b)
+        a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa)
+        b = (a - _sqrt(2.0 * a))/(2.0 * kappa)
+        r = (1.0 + b * b)/(2.0 * b)
 
-    while 1:
-        u1 = random()
+        while 1:
+            u1 = random()
 
-        z = cos(pi * u1)
-        f = (1.0 + r * z)/(r + z)
-        c = kappa * (r - f)
+            z = _cos(_pi * u1)
+            f = (1.0 + r * z)/(r + z)
+            c = kappa * (r - f)
 
-        u2 = random()
+            u2 = random()
 
-        if not (u2 >= c * (2.0 - c) and u2 > c * exp(1.0 - c)):
-            break
+            if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)):
+                break
 
-    u3 = random()
-    if u3 > 0.5:
-        theta = (mu % TWOPI) + acos(f)
-    else:
-        theta = (mu % TWOPI) - acos(f)
+        u3 = random()
+        if u3 > 0.5:
+            theta = (mu % TWOPI) + _acos(f)
+        else:
+            theta = (mu % TWOPI) - _acos(f)
 
-    return theta
+        return theta
 
 # -------------------- gamma distribution --------------------
 
-LOG4 = log(4.0)
-verify('LOG4', 1.38629436111989)
-
-def gammavariate(alpha, beta):
-    # beta times standard gamma
-    ainv = sqrt(2.0 * alpha - 1.0)
-    return beta * stdgamma(alpha, ainv, alpha - LOG4, alpha + ainv)
-
-SG_MAGICCONST = 1.0 + log(4.5)
-verify('SG_MAGICCONST', 2.50407739677627)
-
-def stdgamma(alpha, ainv, bbb, ccc):
-    # ainv = sqrt(2 * alpha - 1)
-    # bbb = alpha - log(4)
-    # ccc = alpha + ainv
-
-    if alpha <= 0.0:
-        raise ValueError, 'stdgamma: alpha must be > 0.0'
-
-    if alpha > 1.0:
-
-        # Uses R.C.H. Cheng, "The generation of Gamma
-        # variables with non-integral shape parameters",
-        # Applied Statistics, (1977), 26, No. 1, p71-74
-
-        while 1:
-            u1 = random()
-            u2 = random()
-            v = log(u1/(1.0-u1))/ainv
-            x = alpha*exp(v)
-            z = u1*u1*u2
-            r = bbb+ccc*v-x
-            if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= log(z):
-                return x
-
-    elif alpha == 1.0:
-        # expovariate(1)
-        u = random()
-        while u <= 1e-7:
+    def gammavariate(self, alpha, beta):
+        # beta times standard gamma
+        ainv = _sqrt(2.0 * alpha - 1.0)
+        return beta * self.stdgamma(alpha, ainv, alpha - LOG4, alpha + ainv)
+
+    def stdgamma(self, alpha, ainv, bbb, ccc):
+        # ainv = sqrt(2 * alpha - 1)
+        # bbb = alpha - log(4)
+        # ccc = alpha + ainv
+
+        random = self.random
+        if alpha <= 0.0:
+            raise ValueError, 'stdgamma: alpha must be > 0.0'
+
+        if alpha > 1.0:
+
+            # Uses R.C.H. Cheng, "The generation of Gamma
+            # variables with non-integral shape parameters",
+            # Applied Statistics, (1977), 26, No. 1, p71-74
+
+            while 1:
+                u1 = random()
+                u2 = random()
+                v = _log(u1/(1.0-u1))/ainv
+                x = alpha*_exp(v)
+                z = u1*u1*u2
+                r = bbb+ccc*v-x
+                if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z):
+                    return x
+
+        elif alpha == 1.0:
+            # expovariate(1)
             u = random()
-        return -log(u)
-
-    else:   # alpha is between 0 and 1 (exclusive)
-
-        # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
-
-        while 1:
-            u = random()
-            b = (e + alpha)/e
-            p = b*u
-            if p <= 1.0:
-                x = pow(p, 1.0/alpha)
-            else:
-                # p > 1
-                x = -log((b-p)/alpha)
-            u1 = random()
-            if not (((p <= 1.0) and (u1 > exp(-x))) or
-                      ((p > 1)  and  (u1 > pow(x, alpha - 1.0)))):
-                break
-        return x
+            while u <= 1e-7:
+                u = random()
+            return -_log(u)
+
+        else:   # alpha is between 0 and 1 (exclusive)
+
+            # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
+
+            while 1:
+                u = random()
+                b = (_e + alpha)/_e
+                p = b*u
+                if p <= 1.0:
+                    x = pow(p, 1.0/alpha)
+                else:
+                    # p > 1
+                    x = -_log((b-p)/alpha)
+                u1 = random()
+                if not (((p <= 1.0) and (u1 > _exp(-x))) or
+                          ((p > 1)  and  (u1 > pow(x, alpha - 1.0)))):
+                    break
+            return x
 
 
 # -------------------- Gauss (faster alternative) --------------------
 
-gauss_next = None
-def gauss(mu, sigma):
-
-    # When x and y are two variables from [0, 1), uniformly
-    # distributed, then
-    #
-    #    cos(2*pi*x)*sqrt(-2*log(1-y))
-    #    sin(2*pi*x)*sqrt(-2*log(1-y))
-    #
-    # are two *independent* variables with normal distribution
-    # (mu = 0, sigma = 1).
-    # (Lambert Meertens)
-    # (corrected version; bug discovered by Mike Miller, fixed by LM)
-
-    # Multithreading note: When two threads call this function
-    # simultaneously, it is possible that they will receive the
-    # same return value.  The window is very small though.  To
-    # avoid this, you have to use a lock around all calls.  (I
-    # didn't want to slow this down in the serial case by using a
-    # lock here.)
-
-    global gauss_next
-
-    z = gauss_next
-    gauss_next = None
-    if z is None:
-        x2pi = random() * TWOPI
-        g2rad = sqrt(-2.0 * log(1.0 - random()))
-        z = cos(x2pi) * g2rad
-        gauss_next = sin(x2pi) * g2rad
-
-    return mu + z*sigma
+    def gauss(self, mu, sigma):
+
+        # When x and y are two variables from [0, 1), uniformly
+        # distributed, then
+        #
+        #    cos(2*pi*x)*sqrt(-2*log(1-y))
+        #    sin(2*pi*x)*sqrt(-2*log(1-y))
+        #
+        # are two *independent* variables with normal distribution
+        # (mu = 0, sigma = 1).
+        # (Lambert Meertens)
+        # (corrected version; bug discovered by Mike Miller, fixed by LM)
+
+        # Multithreading note: When two threads call this function
+        # simultaneously, it is possible that they will receive the
+        # same return value.  The window is very small though.  To
+        # avoid this, you have to use a lock around all calls.  (I
+        # didn't want to slow this down in the serial case by using a
+        # lock here.)
+
+        random = self.random
+        z = self.gauss_next
+        self.gauss_next = None
+        if z is None:
+            x2pi = random() * TWOPI
+            g2rad = _sqrt(-2.0 * _log(1.0 - random()))
+            z = _cos(x2pi) * g2rad
+            self.gauss_next = _sin(x2pi) * g2rad
+
+        return mu + z*sigma
 
 # -------------------- beta --------------------
 
-def betavariate(alpha, beta):
+    def betavariate(self, alpha, beta):
 
-    # Discrete Event Simulation in C, pp 87-88.
+        # Discrete Event Simulation in C, pp 87-88.
 
-    y = expovariate(alpha)
-    z = expovariate(1.0/beta)
-    return z/(y+z)
+        y = self.expovariate(alpha)
+        z = self.expovariate(1.0/beta)
+        return z/(y+z)
 
 # -------------------- Pareto --------------------
 
-def paretovariate(alpha):
-    # Jain, pg. 495
+    def paretovariate(self, alpha):
+        # Jain, pg. 495
 
-    u = random()
-    return 1.0 / pow(u, 1.0/alpha)
+        u = self.random()
+        return 1.0 / pow(u, 1.0/alpha)
 
 # -------------------- Weibull --------------------
 
-def weibullvariate(alpha, beta):
-    # Jain, pg. 499; bug fix courtesy Bill Arms
+    def weibullvariate(self, alpha, beta):
+        # Jain, pg. 499; bug fix courtesy Bill Arms
 
-    u = random()
-    return alpha * pow(-log(u), 1.0/beta)
-
-# -------------------- shuffle --------------------
-# Not quite a random distribution, but a standard algorithm.
-# This implementation due to Tim Peters.
-
-def shuffle(x, random=random, int=int):
-    """x, random=random.random -> shuffle list x in place; return None.
-
-    Optional arg random is a 0-argument function returning a random
-    float in [0.0, 1.0); by default, the standard random.random.
-
-    Note that for even rather small len(x), the total number of
-    permutations of x is larger than the period of most random number
-    generators; this implies that "most" permutations of a long
-    sequence can never be generated.
-    """
-
-    for i in xrange(len(x)-1, 0, -1):
-    # pick an element in x[:i+1] with which to exchange x[i]
-        j = int(random() * (i+1))
-        x[i], x[j] = x[j], x[i]
+        u = self.random()
+        return alpha * pow(-_log(u), 1.0/beta)
 
 # -------------------- test program --------------------
 
-def test(N = 200):
-    print 'TWOPI         =', TWOPI
-    print 'LOG4          =', LOG4
-    print 'NV_MAGICCONST =', NV_MAGICCONST
-    print 'SG_MAGICCONST =', SG_MAGICCONST
-    test_generator(N, 'random()')
-    test_generator(N, 'normalvariate(0.0, 1.0)')
-    test_generator(N, 'lognormvariate(0.0, 1.0)')
-    test_generator(N, 'cunifvariate(0.0, 1.0)')
-    test_generator(N, 'expovariate(1.0)')
-    test_generator(N, 'vonmisesvariate(0.0, 1.0)')
-    test_generator(N, 'gammavariate(0.5, 1.0)')
-    test_generator(N, 'gammavariate(0.9, 1.0)')
-    test_generator(N, 'gammavariate(1.0, 1.0)')
-    test_generator(N, 'gammavariate(2.0, 1.0)')
-    test_generator(N, 'gammavariate(20.0, 1.0)')
-    test_generator(N, 'gammavariate(200.0, 1.0)')
-    test_generator(N, 'gauss(0.0, 1.0)')
-    test_generator(N, 'betavariate(3.0, 3.0)')
-    test_generator(N, 'paretovariate(1.0)')
-    test_generator(N, 'weibullvariate(1.0, 1.0)')
-
-def test_generator(n, funccall):
+def _test_generator(n, funccall):
     import time
     print n, 'times', funccall
     code = compile(funccall, funccall, 'eval')
@@ -356,9 +467,54 @@ def test_generator(n, funccall):
     t1 = time.time()
     print round(t1-t0, 3), 'sec,',
     avg = sum/n
-    stddev = sqrt(sqsum/n - avg*avg)
+    stddev = _sqrt(sqsum/n - avg*avg)
     print 'avg %g, stddev %g, min %g, max %g' % \
               (avg, stddev, smallest, largest)
 
+def _test(N=200):
+    print 'TWOPI         =', TWOPI
+    print 'LOG4          =', LOG4
+    print 'NV_MAGICCONST =', NV_MAGICCONST
+    print 'SG_MAGICCONST =', SG_MAGICCONST
+    _test_generator(N, 'random()')
+    _test_generator(N, 'normalvariate(0.0, 1.0)')
+    _test_generator(N, 'lognormvariate(0.0, 1.0)')
+    _test_generator(N, 'cunifvariate(0.0, 1.0)')
+    _test_generator(N, 'expovariate(1.0)')
+    _test_generator(N, 'vonmisesvariate(0.0, 1.0)')
+    _test_generator(N, 'gammavariate(0.5, 1.0)')
+    _test_generator(N, 'gammavariate(0.9, 1.0)')
+    _test_generator(N, 'gammavariate(1.0, 1.0)')
+    _test_generator(N, 'gammavariate(2.0, 1.0)')
+    _test_generator(N, 'gammavariate(20.0, 1.0)')
+    _test_generator(N, 'gammavariate(200.0, 1.0)')
+    _test_generator(N, 'gauss(0.0, 1.0)')
+    _test_generator(N, 'betavariate(3.0, 3.0)')
+    _test_generator(N, 'paretovariate(1.0)')
+    _test_generator(N, 'weibullvariate(1.0, 1.0)')
+
+# Initialize from current time.
+_inst = Random()
+seed = _inst.seed
+random = _inst.random
+uniform = _inst.uniform
+randint = _inst.randint
+choice = _inst.choice
+randrange = _inst.randrange
+shuffle = _inst.shuffle
+normalvariate = _inst.normalvariate
+lognormvariate = _inst.lognormvariate
+cunifvariate = _inst.cunifvariate
+expovariate = _inst.expovariate
+vonmisesvariate = _inst.vonmisesvariate
+gammavariate = _inst.gammavariate
+stdgamma = _inst.stdgamma
+gauss = _inst.gauss
+betavariate = _inst.betavariate
+paretovariate = _inst.paretovariate
+weibullvariate = _inst.weibullvariate
+getstate = _inst.getstate
+setstate = _inst.setstate
+
 if __name__ == '__main__':
-    test()
+    _test()
diff --git a/Misc/NEWS b/Misc/NEWS
index e48add5..058d55c 100644
--- a/Misc/NEWS
+++ b/Misc/NEWS
@@ -1,3 +1,16 @@
+What's New in Python 2.1 alpha 2?
+=================================
+Core language, builtins, and interpreter
+
+
+Standard library
+
+- random.py is now self-contained, and offers all the functionality of
+  the now-deprecated whrandom.py.  See the docs for details.  random.py
+  also supports new functions getstate() and setstate(), for saving
+  and restoring the internal state of all the generators.
+
+
 What's New in Python 2.1 alpha 1?
 =================================