from random import random, shuffle
from math import log, sqrt, exp

LOG4 = log(4.0)
SG_MAGICCONST = 1.0 + log(4.5)

def gammavariate(alpha, beta):
    if alpha > 1.0:
        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):
                return x * beta

    elif alpha == 1.0:
        return -log(1.0 - random()) * beta

    else:
        while True:
            u = random()
            b = (_e + alpha)/_e
            p = b*u
            if p <= 1.0:
                x = p ** (1.0/alpha)
            else:
                x = -log((b-p)/alpha)
            u1 = random()
            if p > 1.0:
                if u1 <= x ** (alpha - 1.0):
                    break
            elif u1 <= exp(-x):
                break
        return x * beta

def betavariate(alpha, beta):
    y = gammavariate(alpha, 1.0)
    if y == 0:
        return 0.0
    else:
        return y / (y + gammavariate(beta, 1.0))

def binomialvariate(n, p):
    if n == 0:
        return 0
    a, b = (1 + n)//2, 1 + n//2
    x = betavariate(a, b)
    if x < p:
        return a + binomialvariate(b - 1, (p - x)/(1 - x))
    else:
        return binomialvariate(a - 1, p/x)

def choices(population, weights, *, k=1):
    r = 1.0
    n = k
    selections = []
    for elem, p in zip(population, weights):
        v = binomialvariate(n, p / r)
        selections += [elem] * v
        n -= v
        r -= p
    shuffle(selections)
    return selections

if __name__ == '__main__':

    from collections import Counter

    print(sorted(Counter(choices('ABCD', (0.10, 0.30, 0.20, 0.40), k=100_000)).items()))