from scipy.misc import comb


def hcp(i, j, k, l):
    return (4 * i) + (3 * j) + (2 * k) + l


def var(i, j, k, l):
    return (hcp(i, j, k, l) - 10)**2


def expectation(func):
    mean = 0
    for i in range(5):
        I = comb(4, i, exact=True)
        for j in range(5):
            J = comb(4, j, exact=True)
            for k in range(5):
                K = comb(4, k, exact=True)
                for l in range(5):
                    L = comb(4, l, exact=True)
                    total = i + j + k + l
                    if total <= 13:
                        p = I * J * K * L * comb(36, 13-total, exact=True) / float(comb(52, 13, exact=True))
                        mean += func(i, j, k, l) * p
    return mean


def main():
    print "Expectation is " + str(expectation(hcp))
    print "Variance is " + str(expectation(var))

if __name__ == '__main__':
    main()