import random

def roll_dice(n):
    """Throw a dice n times, count number of sixes"""
    M = 0
    for i in range(n):
        dice = random.randint(1,6)
        if dice == 6:
            M += 1
    return M

def all_sixes(n,N=1000):
    M = 0
    for i in range(N):
        result = roll_dice(n)
        if result == n:
            M += 1
    return float(M)/N


def next_six(n,N=1000):
    M1 = 0
    M2 = 0

    for i in range(N):
        if roll_dice(n) == n:
            M1 += 1
            result = roll_dice(1)
            if result == 1:
                M2 += 1

    return float(M2)/M1



print("Estimated probability of throwing 4 sixes: %g, analytical probability: %g" %(all_sixes(4,10000), 1.0/6**4))
print("Probability of next six ", next_six(3))

"""
NOTE: Although this is a correctly implemented Monte Carlo simulation for this problem,
the question is not very well posed. If there are many dice, the probability of
getting all sixes is very low, and therefore the number of times we throw the
last dice in next_six is very low, and the results are variable and inaccurate. We
can solve the problem by choosing N very large, but then the problem takes a long
time to solve."""