Solution to a problem posed by Deomani Sharma Ramsurrun

gistfile1.txt

# Solution to a problem posted by Deomani Sharma Ramsurrun
#       David Radcliffe  11 Dec 2013
# This code is in the public domain
###########################################################

import fractions

# Generate all non-negative integer solutions to the system
#   ab + ac + ad = 6
#   ab + bc + bd = 9
#   ac + bc + cd = 13 
#   ad + bd + cd = 14

solutions = [(ab, ac, 6-ab-ac, 7-ab-ac, ac+2, ab+6)
             for ab in range(7)
             for ac in range(7-ab)]

# Multinomial coefficient
#   multinomial((a,b,c,d)) = (a+b+c+d)! / (a! * b! * c! * d!)

def multinomial(vec):
    v = list(vec)
    i = 0
    N = sum(v)
    result = fractions.Fraction(1)
    while i < len(vec) - 1:
        if v[i] > 0:
            result = result*N/v[i]
            N -= 1
            v[i] -= 1
        else:
            i += 1
    return int(result)

# Count the number of ways that each solution can be generated

counts = map(multinomial, solutions)
numer = [sum(s[i]*c for s,c in zip(solutions, counts)) for i in range(6)]
denom = sum(numer)
labels = ("AB", "AC", "AD", "BC", "BD", "CD")
for i in range(6):
    frac = fractions.Fraction(numer[i], denom)
    print "P(%s) = %s = %f " % (labels[i], frac, frac)

# Expected output:
#
# P(AB) = 68378/1129541 = 0.060536 
# P(AC) = 351136/3388623 = 0.103622 
# P(AD) = 58844/484089 = 0.121556 
# P(BC) = 573271/3388623 = 0.169175 
# P(BD) = 96266/484089 = 0.198860 
# P(CD) = 55872/161363 = 0.346250