def binprob(n, phi, integer=True):
    """
    Computes all binomial terms for n given the Bernoulli probability phi. 
    Returns the relative frequencies AND the cumulative frequencies (two lists). 
    """

    assert is_posinteger(n),  \
                      "first argument to binProb must be a positive integer!"
    assert 0.0 <= phi <= 1.0, \
                    "Bernoulli probability must be in [0.0, 1.0] in binProb!"

    farray = []
    np1    = n + 1
    for k in range(0, np1):
        x = fbincoeff(n, k, integer)
        farray.append(x)

    parray = []
    carray = []

    y = 0.0
    q = 1.0 - phi
    for k in range(0, np1):
        x  = farray[k] * phi**k * q**(n-k)
        x  = kept_within(0.0, x, 1.0)
        y  = y + x
        y  = kept_within(0.0, y, 1.0)
        parray.append(x)
        carray.append(y)

    return parray, carray

def ipareto_zero(prob, lam, xm=1.0):
    """
    The inverse of the Pareto distribution with the support shifted to [0, inf):
    f = lam * xm**lam / (x+xm)**(lam+1)
    F = 1 - [xm/(x+xm)]**lam
    x in [0, inf)
    lam > 0
    For lam < 1 all moments are infinite
    For lam < 2 all moments are infinite except for the mean
    """

    _assertprob(prob, 'ipareto_zero')
    assert lam > 0.0, "shape parameter lambda in ipareto_zero must be positive!"
    textxm1 = "left support limit parameter xm of unshifted in ipareto_zero"
    textxm2 = "distribution must not be negative in ipareto_zero!"
    assert xm  >= 0.0, textxm1 + textxm2

    q  =  1.0 - prob

    if q == 0.0: return float('inf')

    x  =  xm * (safepow(q, -1.0/lam) - 1.0)

    x  =  kept_within(0.0, x)

    return x

def ipareto(prob, lam, xm=1.0):
    """
    The inverse of the Pareto distribution: 
    f = lam * xm**lam / x**(lam+1) 
    F = 1 - (xm/x)**lam
    x in [xm, inf)
    lam > 0
    For lam < 1 all moments are infinite
    For lam < 2 all moments are infinite except for the mean
    """

    _assertprob(prob, 'ipareto')
    assert lam >  0.0, "shape parameter lambda in ipareto must be positive!"
    assert xm  >= 0.0, \
          "left support limit parameter xm must not be negative in ipareto!"

    q  =  1.0 - prob

    if q == 0.0: return float('inf')
    
    x  =  xm * safepow(q, -1.0/lam)

    x  =  kept_within(xm, x)

    return x

    def rtri_unif_tri(self, a, b, c, d):
        """
        Triangular-uniform-triangular distribution with support on [a, d] and 
        with breakpoints in b and c
                      ------
        pdf:        /        \
                   /           \
            ------               -------                                         
        """

        # Input check -----------------------
        assert a <= b and b <= c and c <= d, \
                                "break points scrambled in rtri_unif_tri!"
        # -----------------------------------


        if d == a: return a


        dcba   =  d + c - b - a
        h      =  2.0 / dcba
        first  =  0.5 * h * (b-a)
        p      =  self.runif01()
        poh    =  0.5 * p * dcba

        if p <= first:
            x  =  sqrt(2.0*(b-a)*poh) + a
        elif first < p <= first + h*(c-b):
            x  =  (c-b)*(poh-0.5*(b-a)) + b
        else:
            x  =  d - sqrt((d-c)*dcba*(1.0-p))

        x  = kept_within(a, x, d)

        return x

    def rtriang(self, left, mode, right):
        """
        Generator of triangularly distributed random numbers on [left, right] 
        with the peak of the pdf at mode. 
        """

        assert left <= mode and mode <= right, \
                                  "mode out of support range in rtriang!"

        p  =  self.runif01()

        span    =  right - left
        spanlo  =  mode  - left
        spanhi  =  right - mode
        #height  =  2.0 / span
        #surf1   =  0.5 * spanlo * height
        #surf1   =  spanlo/float(span)

        #if p <= surf1:
        if p <= spanlo/float(span):
            #x  =  sqrt(2.0*spanlo*p/height)
            x  =  sqrt(spanlo*span*p)
        else:
            #x  =  span - sqrt(2.0*spanhi*(1.0-p)/height)
            x  =  span - sqrt(spanhi*span*(1.0-p))
        x += left

        x  = kept_within(left, x, right)

        return x

def icoxian2(prob, means, probs):
    """
    The Coxian phased distribution, which is based on the exponential.
    probs is a list of probabilities for GOING ON TO THE NEXT PHASE rather 
    than reaching the absorbing state prematurely. The number of means must 
    (of course) be one more than the number of probabilities! 
    
    NB No two means[k] must be equal - if equal means are is desired, use 
    icoxian instead (slower, however). 
    """

    _assertprob(prob, 'icoxian2')
    # Everything else will be checked in dcoxian2 and ccoxian2

   # ------------------------------------
    def _fifi2fid(x):
        x      = kept_within(0.0, x)
        cdf    = ccoxian2(means, probs, x)
        pdf    = dcoxian2(means, probs, x)
        fi     = cdf - prob
        if pdf <= 0.0:
            if fi == 0.0: fi2fid = 1.0
            else:         fi2fid = MAXFLOAT
        else:
            fi2fid = fi/pdf
        return fi, fi2fid
   # ------------------------------------

    x = znewton(_fifi2fid, 0.0, 'icoxian2', tolf=SQRTMACHEPS)

    x = kept_within(0.0, x)
    
    return x

def itriang(prob, left, mode, right):
    """
    The inverse of a triangular distribution between left and right 
    having its peak at x = mode 
    """

    _assertprob(prob, 'itriang')
    # ---
    assert left <= mode and mode <= right, \
                                  "mode out of support range in itriang!"

    span    =  right - left
    spanlo  =  mode  - left
    spanhi  =  right - mode
    #height  =  2.0 / span
    #surf1   =  0.5 * spanlo * height
    #surf1   =  spanlo/float(span)

    #if prob <= surf1:
    if prob <= spanlo/float(span):
        #x  =  sqrt(2.0*spanlo*prob/height)
        x  =  sqrt(spanlo*span*prob)
    else:
        #x  =  span - sqrt(2.0*spanhi*(1.0-prob)/height)
        x  =  span - sqrt(spanhi*span*(1.0-prob))
    x += left

    x  = kept_within(left, x, right)

    return x

def ibinomial(prob, n, phi, normconst=10.0):
    """
    The binomial distribution: p(N=k) = bincoeff * phi**k * (1-phi)**(n-k), 
    n >= 1, k = 0, 1,...., n 
    """

    # Input check -----------
    _assertprob(prob, 'ibinomial')
    assert is_posinteger(n),        "n must be a positive integer in ibinomial!"
    assert 0.0 <= phi and phi <= 1.0, \
                          "success frequency out of support range in ibinomial!"
    assert normconst >= 10.0, \
           "parameter limit for normal approx. in ibinomial must not be < 10.0!"
    # -----------------------

    onemphi = 1.0 - phi

    if phi < 0.5: w = normconst * onemphi / phi
    else:         w = normconst * phi / onemphi

    if n > w:
        k = int(round(inormal(prob, n*phi, sqrt(n*phi*onemphi))))

    else:
        k   = 0
        cdf = binProb(n, phi)[1]
        while True:
            if cdf[k] <= prob:  k = k + 1
            else:               break

    k = kept_within(0, k, n)

    return k

    def rbinomial(self, n, phi):
        """
        The binomial distribution: p(N=k) = bincoeff * phi**k * (1-phi)**(n-k); 
        n >= 1;  k = 0, 1,...., n  where phi is the frequency or "Bernoulli 
        probability".
        
        Algorithm taken from ORNL-RSIC-38, Vol II (1973). 
        """

        assert is_posinteger(n), \
                          "n must be a positive integer in rbinomial!"
        assert 0.0 < phi and phi < 1.0, \
                   "frequency parameter is out of range in rbinomial!"

        normconst = 10.0
        onemphi   = 1.0 - phi
        if phi < 0.5: w = int(round(normconst * onemphi / phi))
        else:         w = int(round(normconst * phi / onemphi))

        if n > w:
            #-------------------------------------------------------
            k = int(round(self.rnormal(n*phi, sqrt(n*phi*onemphi))))

        else:
            #-------------------------------------------------------
            if phi < 0.25:
                k   = -1
                m   =  0
                phi = - safelog(onemphi)
                while m < n:
                    r  = self.rexpo(1.0)
                    j  = 1 + int(r/phi)
                    m += j
                    k += 1
                if m == n:
                    k += 1

            elif phi > 0.75:
                k   =  n + 1
                m   =  0
                phi = - safelog(phi)
                while m < n:
                    r  = self.rexpo(1.0)
                    j  = 1 + int(r/phi)
                    m += j
                    k -= 1
                if m == n:
                    k -= 1

            else: # if 0.25 <= phi and phi <= 0.75:
                k = 0
                m = 0
                while m < n:
                    r  = self.runif01()
                    if r < phi: k += 1
                    m += 1

        k = kept_within(0, k, n)

        return k

def cexppower(loc, scale, alpha, x, lngam1oalpha=False, tolf=FOURMACHEPS, itmax=128):
    """
    The exponential power distribution 
    f  =  (a/s) * exp(-abs([x-l]/s)**a) / [2*gamma(1/a)]
    F  =  1/2 * [1 + sgn(x-l) * Fgamma(1/a, abs([x-l]/s)**a)],   x in R
    s, a > 0
    where Fgamma is the gamma distribution cdf.

    NB It is possible to gain efficiency by providing the value of the 
    natural logarithm of the complete gamma function ln(gamma(1.0/alpha)) 
    as a pre-computed input (may be computed using numlib.specfunc.lngamma) 
    instead of the default 'False'.

    tolf and itmax are the numerical control parameters of cgamma.
    """

    assert scale > 0.0, "scale parameter must be a positive float in cexppower!"
    assert alpha > 0.0, "shape parameter alpha must be a positive float in cexppower!"

    if alpha == 1.0:
        return claplace(loc, scale, x)

    ainv = 1.0 / alpha
    xml = x - loc

    if not lngam1oalpha:
        lng1oa = lngamma(ainv)
    else:
        lng1oa = lngam1oalpha
    cg = cgamma(ainv, 1.0, abs(xml / scale) ** alpha, lng1oa, tolf, itmax)
    cdf = 0.5 * (fsign(xml) * cg + 1.0)

    cdf = kept_within(0.0, cdf, 1.0)

    return cdf

def ctriang(left, mode, right, x):
    """
    The cdf of the triangular distribution with support 
    on [left, right] and with mode 'mode'. 
    """

    # Input check -----------------------
    assert right > left, "support range must be positive in ctriang!"
    assert left <= mode and mode <= right, "mode must be within support range in ctriang!"
    assert left <= x and x <= right, "variate must be within support range in ctriang!"
    # -----------------------------------

    spant = right - left
    spanl = mode - left
    spanr = right - mode

    if spanr == 0.0:
        cdf = (x - left) ** 2 / float((spant * spanl))

    elif spanl == 0.0:
        cdf = 1.0 - (right - x) ** 2 / float((spant * spanr))

    elif x <= mode:
        cdf = (x - left) ** 2 / float((spant * spanl))

    else:
        cdf = 1.0 - (right - x) ** 2 / float((spant * spanr))

    cdf = kept_within(0.0, cdf, 1.0)

    return cdf

def cerlang(nshape, phasemean, x):
    """
    The cdf of the Erlang distribution.
    Represents the sum of nshape exponentially distributed random variables, 
    all having "phasemean" as mean
    """

    if nshape == 1:
        cdf = cexpo(phasemean, x)

    else:
        assert is_posinteger(nshape), "shape parameter must be a positive integer in cerlang!"
        assert phasemean >= 0.0, "phase mean must not be negative in cerlang!"
        assert x >= 0.0, "variate must not be negative in cerlang!"
        y = x / float(phasemean)
        cdf = 1.0
        term = 1.0
        cdf = term
        for k in range(1, nshape):
            term = term * y / k
            cdf = cdf + term

        cdf = 1.0 - exp(-y) * cdf

        cdf = kept_within(0.0, cdf, 1.0)

    return cdf

def clevy(scale, x):
    """
    The cdf of the Levy distribution (stable distribution with 
    alpha = 1/2 and beta = 1, aka the Cournot distribution). 
    This is actually the right-skewed Levy!
    f = sqrt(s/2pi) * (1/x)**(3/2) * exp(-s/2x)
    F = erfc(sqrt(s/2x))
    
    s >= 0.0, x >= 0
    """

    assert scale >= 0.0, "scale must not be negative in clevy!"
    assert x >= 0.0, "variate must not be negative in clevy!"

    # The cdf of the Levy can be handled since it is an "incomplete gamma
    # function", but it seems to be more simple than that:

    try:
        cdf = erfc1(sqrt(0.5 * scale / x))
    except (OverflowError, ZeroDivisionError):
        return 0.0

    cdf = kept_within(0.0, cdf, 1.0)

    return cdf

def _stable_sym_tail(alpha, x):
    """
    An asymptotic expression for the tail.
    """

    # calpha = exp(lngamma(alpha)) * sin(PIHALF*alpha) / PI
    calpha = PIINV * exp(lngamma(alpha)) * sin(PIHALF * alpha)

    try:
        cdf = calpha / x ** alpha

    except ZeroDivisionError:
        cdf = log(calpha) - alpha * log(x)
        try:
            cdf = exp(cdf)
        except OverflowError:
            cdf = 0.0

    except OverflowError:
        cdf = log(calpha) - alpha * log(x)
        try:
            cdf = exp(cdf)
        except OverflowError:
            cdf = 0.0

    cdf = 1.0 - cdf

    cdf = kept_within(0.5, cdf, 1.0)
    return cdf

def cextreme_I(type, mu, scale, x):
    """
    Extreme value distribution type I (aka the Gumbel distribution or 
    Gumbel distribution type I):
    F = exp{-exp[-(x-mu)/scale]}       (max variant)
    f = exp[-(x-mu)/scale] * exp{-exp[-(x-mu)/scale]} / scale
    F = 1 - exp{-exp[+(x-mu)/scale]}   (min variant)
    f = exp[+(x-mu)/scale] * exp{-exp[+(x-mu)/scale]} / scale

    type must be 'max' or 'min'
    scale must be > 0.0
    """

    assert scale > 0.0, "scale must be positive in cextreme_I!"

    if type == "max":
        cdf = exp(-exp(-(x - mu) / float(scale)))
    elif type == "min":
        cdf = 1.0 - exp(-exp((x - mu) / float(scale)))
    else:
        raise Error("type must be either 'max' or 'min' in cextreme_I!")

    cdf = kept_within(0.0, cdf, 1.0)

    return cdf

def ibeta(prob, a, b, x1=0.0, x2=1.0, betaab=False):
    """
    The beta distribution:
    f = x**(a-1) * (1-x)**(b-1) / beta(a, b)
    a, b >= 0; 0 <= x <= 1
    F is the integral = the incomplete beta or the incomplete beta ratio 
    function depending on how the incomplete beta function is defined.
    
    x2 >= x1 !!!!
    
    NB It is possible to provide the value of the complete beta 
    function beta(a, b) as a pre-computed input (may be computed 
    using numlib.specfunc.beta) instead of the default "False", 
    a feature that will make ibeta 30 % faster!
    """

    # Everything will be checked in cbeta

    if a == 1.0 and b == 1.0: return iunifab(prob, x1, x2)

   # -----------------------------------------------------------
    def _fi(x):
        return cbeta(a, b, x1, x2, x, betaab) - prob
   # -----------------------------------------------------------

    x = zbrent(_fi, x1, x2, 'ibeta', tolf=SQRTMACHEPS)

    x = kept_within(x1, x, x2)
    
    return x

    def rerlang(self, nshape, phasemean, xmax=float('inf')):
        """
        Generator of Erlang-distributed random variates.
        Represents the sum of nshape exponentially distributed random variables, 
        each having the same mean value = phasemean. For nshape = 1 it works as 
        a generator of exponentially distributed random numbers.
        """

        assert is_posinteger(nshape), \
                       "shape parameter must be a positive integer in rerlang!"
        assert phasemean >= 0.0,   "phasemean must not be negative in rerlang!"
        assert xmax >= 0.0,      "variate max must be non-negative in rerlang!"


        if nshape < GeneralRandomStream.__ERLANG2GAMMA:
            while True:
                x  =  1.0
                for k in range(0, nshape):
                    x *= self.runif01() # Might turn out to be zero...
                x  = - phasemean * safelog(x)
                if x <= xmax: break

        else:   # Gamma is OK
            while True:
                x  =  phasemean * self.rgamma(float(nshape), 1.0)
                if x <= xmax: break

        x  =  kept_within(0.0, x)

        return x

def iNexpo2(prob, means):
    """
    A distribution of a sum of exponential random variables.
    
    NB No two means are allowed to be equal!!!! 
    """

    _assertprob(prob, 'iNexpo2')
    # Everything else will be checked in iexpo, cNexpo2 and dNexpo2

    if len(means) == 1:
        x = iexpo(prob, means[0])

    else:
       # ------------------------------------
        def _fifi2fid(x):
            x      = kept_within(0.0, x)
            cdf    = cNexpo2(means, x)
            pdf    = dNexpo2(means, x)
            fi     = cdf - prob
            if pdf <= 0.0:
                if fi == 0.0: fi2fid = 1.0
                else:         fi2fid = MAXFLOAT
            else:
                fi2fid = fi/pdf
            return fi, fi2fid
       # ------------------------------------

        x = znewton(_fifi2fid, max(means), 'iNexpo2', tolf=SQRTMACHEPS)

        x = kept_within(0.0, x)
    
    return x

    def antithet_sample(self, nparams):
        """
        Generates a matrix having two rows, the first row being a list of 
        uniformly distributed random numbers p in [0.0, 1.0], each row 
        containing nparams elements. The second row contains the corresponding 
        antithetic sample with the complements 1-p. 
        """

        rstream = self.rstream

        antimatrix = Matrix()  # antimatrix belongs to the Matrix class
        for k in range(0, nparams):
            pvector = array('d', [])
            p1  =  rstream.runif01()
            pvector.append(p1)
            dum =  rstream.runif01()  # For synchronization only - never used
            p2  =  1.0 - p1
            p2 = kept_within(0.0, p2, 1.0) # Probabilities must be in [0.0, 1.0]
            pvector.append(p2)
            antimatrix.append(pvector)

        # Matrix must be transposed in order for each sample to occupy one row.
        # Sample vector k is in antimatrix[k], where k is 0 or 1
        antimatrix.transpose()

        return antimatrix

def ierlang(prob, nshape, phasemean=1.0):
    """
    Represents the sum of nshape exponentially distributed random variables, 
    each having the same mean value = phasemean 
    """

    _assertprob(prob, 'ierlang')
    # Everything else will be checked in iexpo, cerlang and derlang

    if nshape == 1:
        x = iexpo(prob, phasemean)

    else:
       # ------------------------------------
        def _fifi2fid(x):
            x      = kept_within(0.0, x)
            cdf    = cerlang(nshape, phasemean, x)
            pdf    = derlang(nshape, phasemean, x)
            fi     = cdf - prob
            if pdf <= 0.0:
                if fi == 0.0: fi2fid = 1.0
                else:         fi2fid = MAXFLOAT
            else:
                fi2fid = fi/pdf
            return fi, fi2fid
       # ------------------------------------

        x = znewton(_fifi2fid, (nshape-1.0)*phasemean, 'ierlang', \
                                                     tolf=SQRTMACHEPS)

        x = kept_within(0.0, x)
    
    return x