from random import _inst
from math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil
from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
from itertools import accumulate as _accumulate
from operator import add as _add
from bisect import bisect as _bisect

NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
TWOPI = 2.0*_pi
LOG4 = _log(4.0)
SG_MAGICCONST = 1.0 + _log(4.5)
BPF = 53        # Number of bits in a float
RECIP_BPF = 2**-BPF

def random(*, random=_inst):
    return iter(random.random, None)

## -------------------- integer distributions  -------------------

def randrange(start, stop=None, step=1, *, random=_inst):
    """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' argument.
    """

    # This code is a bit messy to make it fast for the
    # common case while still doing adequate error checking.
    randbelow = random._randbelow
    istart = int(start)
    if istart != start:
        raise ValueError("non-integer arg 1 for randrange()")
    if stop is None:
        if istart > 0:
            while True:
                yield randbelow(istart)
        raise ValueError("empty range for randrange()")

    # stop argument supplied.
    istop = int(stop)
    if istop != stop:
        raise ValueError("non-integer stop for randrange()")
    width = istop - istart
    if step == 1 and width > 0:
        while True:
            yield istart + randbelow(width)
    if step == 1:
        raise ValueError("empty range for randrange() (%d,%d, %d)" % (istart, istop, width))

    # Non-unit step argument supplied.
    istep = int(step)
    if istep != step:
        raise ValueError("non-integer step for randrange()")
    if istep > 0:
        n = (width + istep - 1) // istep
    elif istep < 0:
        n = (width + istep + 1) // istep
    else:
        raise ValueError("zero step for randrange()")

    if n <= 0:
        raise ValueError("empty range for randrange()")

    while True:
        yield istart + istep * randbelow(n)

def randint(a, b, *, random=_inst):
    """Return random integer in range [a, b], including both end points.
    """
    return randrange(a, b + 1, random=random)

## -------------------- sequence methods  -------------------

def choice(seq, weights=None, *, random=_inst):
    """Choose a random element from a non-empty sequence.

    Optional argument weights is a sequence of non-negative numbers used to
    alters the probability that element at index i in the sequence is
    selected.
    """
    if weights is None:
        n = len(seq)
        if not n:
            raise IndexError('Cannot choose from an empty sequence')
        randbelow = random._randbelow
        while True:
            yield seq[randbelow(n)]

    cumulative_dist = list(_accumulate(weights, _add))
    total_sum = cumulative_dist[-1]
    if len(weights) != len(seq):
        raise ValueError("Length of weights must equal length of sequence")
    if any(w < 0 for w in weights):
        raise ValueError("All weights must be non-negative")
    if not total_sum:
        raise ValueError("At least one weight must be greater than zero")

    if isinstance(total_sum, int):
        randbelow = random._randbelow
        while True:
            yield seq[_bisect(cumulative_dist, randbelow(total_sum))]
    else:
        for u in iter(random.random, None):
            yield seq[_bisect(cumulative_dist, u * total_sum)]

## -------------------- real-valued distributions  -------------------

## -------------------- uniform distribution -------------------

def uniform(a, b, *, random=_inst):
    "Get a random number in the range [a, b) or [a, b] depending on rounding."
    d = b - a
    for u in iter(random.random, None):
        yield a + d * u

## -------------------- triangular --------------------

def triangular(low=0.0, high=1.0, mode=None, *, random=_inst):
    """Triangular distribution.

    Continuous distribution bounded by given lower and upper limits,
    and having a given mode value in-between.

    http://en.wikipedia.org/wiki/Triangular_distribution

    """
    d = high - low
    c = 0.5 if mode is None else (mode - low) / d
    c2 = 1.0 - c
    for u in iter(random.random, None):
        if u > c:
            yield high - d * ((1.0 - u) * c2) ** 0.5
        else:
            yield low + d * (u * c) ** 0.5

## -------------------- normal distribution --------------------

def normalvariate(mu, sigma, *, random=_inst):
    """Normal distribution.

    mu is the mean, and sigma is the standard deviation.

    """
    # 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 = random.random
    while True:
        u1 = random()
        u2 = 1.0 - random()
        z = NV_MAGICCONST * (u1 - 0.5) / u2
        zz = z * z / 4.0
        if zz <= -_log(u2):
            yield mu + z * sigma

## -------------------- lognormal distribution --------------------

def lognormvariate(mu, sigma, *, random=_inst):
    """Log normal distribution.

    If you take the natural logarithm of this distribution, you'll get a
    normal distribution with mean mu and standard deviation sigma.
    mu can have any value, and sigma must be greater than zero.

    """
    return map(_exp, normalvariate(mu, sigma, random=random))

## -------------------- exponential distribution --------------------

def expovariate(lambd, *, random=_inst):
    """Exponential distribution.

    lambd is 1.0 divided by the desired mean.  It should be
    nonzero.  (The parameter would be called "lambda", but that is
    a reserved word in Python.)  Returned values range from 0 to
    positive infinity if lambd is positive, and from negative
    infinity to 0 if lambd is negative.

    """
    # lambd: rate lambd = 1/mean
    # ('lambda' is a Python reserved word)

    # we use 1-random() instead of random() to preclude the
    # possibility of taking the log of zero.
    for u in iter(random.random, None):
        yield -_log(1.0 - u) / lambd

## -------------------- von Mises distribution --------------------

def vonmisesvariate(mu, kappa, *, random=_inst):
    """Circular data distribution.

    mu is the mean angle, expressed in radians between 0 and 2*pi, and
    kappa is the concentration parameter, which must be greater than or
    equal to zero.  If kappa is equal to zero, this distribution reduces
    to a uniform random angle over the range 0 to 2*pi.

    """
    # 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.

    # Thanks to Magnus Kessler for a correction to the
    # implementation of step 4.

    if kappa <= 1e-6:
        for u in iter(random.random, None):
            yield TWOPI * u

    s = 0.5 / kappa
    r = s + _sqrt(1.0 + s * s)

    random = random.random
    while True:
        u1 = random()
        z = _cos(_pi * u1)

        d = z / (r + z)
        u2 = random()
        if u2 < 1.0 - d * d or u2 <= (1.0 - d) * _exp(d):
            q = 1.0 / r
            f = (q + z) / (1.0 + q * z)
            u3 = random()
            if u3 > 0.5:
                theta = (mu + _acos(f)) % TWOPI
            else:
                theta = (mu - _acos(f)) % TWOPI

            yield theta

## -------------------- gamma distribution --------------------

def gammavariate(alpha, beta, *, random=_inst):
    """Gamma distribution.  Not the gamma function!

    Conditions on the parameters are alpha > 0 and beta > 0.

    The probability distribution function is:

                x ** (alpha - 1) * math.exp(-x / beta)
        pdf(x) =  --------------------------------------
                    math.gamma(alpha) * beta ** alpha

    """

    # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2

    # Warning: a few older sources define the gamma distribution in terms
    # of alpha > -1.0
    if alpha <= 0.0 or beta <= 0.0:
        raise ValueError('gammavariate: alpha and beta must be > 0.0')

    random = random.random
    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

        ainv = _sqrt(2.0 * alpha - 1.0)
        bbb = alpha - LOG4
        ccc = alpha + ainv

        while True:
            u1 = random()
            if not 1e-7 < u1 < .9999999:
                continue
            u2 = 1.0 - 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):
                yield x * beta

    elif alpha == 1.0:
        # expovariate(1)
        for u in iter(random.random, None):
            yield -_log(1.0 - u) * beta

    else:   # alpha is between 0 and 1 (exclusive)

        # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle

        b = (_e + alpha) / _e
        d = 1.0 / alpha
        d2 = alpha - 1.0
        while True:
            u = random()
            p = b * u
            if p <= 1.0:
                x = p ** d
                if random() <= _exp(-x):
                    yield x * beta
            else:
                x = -_log((b - p) * d)
                if random() <= x ** d2:
                    yield x * beta

## -------------------- Gauss (faster alternative) --------------------

def gauss(mu, sigma, *, random=_inst):
    """Gaussian distribution.

    mu is the mean, and sigma is the standard deviation.  This is
    slightly faster than the normalvariate() function.

    Not thread-safe without a lock around calls.

    """

    # 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 = random.random
    while True:
        x2pi = random() * TWOPI
        g2rad = _sqrt(-2.0 * _log(1.0 - random())) * sigma
        yield mu + _cos(x2pi) * g2rad
        yield mu + _sin(x2pi) * g2rad

## -------------------- beta --------------------
## See
## http://mail.python.org/pipermail/python-bugs-list/2001-January/003752.html
## for Ivan Frohne's insightful analysis of why the original implementation:
##
##    def betavariate(alpha, beta):
##        # Discrete Event Simulation in C, pp 87-88.
##
##        y = self.expovariate(alpha)
##        z = self.expovariate(1.0/beta)
##        return z/(y+z)
##
## was dead wrong, and how it probably got that way.

def betavariate(alpha, beta, *, random=_inst):
    """Beta distribution.

    Conditions on the parameters are alpha > 0 and beta > 0.
    Returned values range between 0 and 1.

    """

    # This version due to Janne Sinkkonen, and matches all the std
    # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
    b = gammavariate(beta, 1.0, random=random)
    for y, x in zip(gammavariate(alpha, 1.0, random=random),
                    gammavariate(beta, 1.0, random=random)):
        if y == 0.0:
            yield 0.0
        else:
            yield y / (y + x)

## -------------------- Pareto --------------------

def paretovariate(alpha, *, random=_inst):
    """Pareto distribution.  alpha is the shape parameter."""
    # Jain, pg. 495

    d = -1.0 / alpha
    for u in iter(random.random, None):
        yield (1.0 - u) ** d

## -------------------- Weibull --------------------

def weibullvariate(alpha, beta, *, random=_inst):
    """Weibull distribution.

    alpha is the scale parameter and beta is the shape parameter.

    """
    # Jain, pg. 499; bug fix courtesy Bill Arms

    d = -1.0 / beta
    for u in iter(random.random, None):
        yield alpha * (-_log(1.0 - u)) ** d