import json, math, random
from collections import deque
import numpy as np

SEED = 418
random.seed(SEED); np.random.seed(SEED)

# Cyclic group Z_m, represented additively. Its 2-D orthogonal representation is a rotation.
def inv(g, m): return (-g) % m
def compose(a, b, m): return (a + b) % m
def rotation(g, m):
    t = 2 * math.pi * g / m
    return np.array([[math.cos(t), -math.sin(t)], [math.sin(t), math.cos(t)]])

def add_edge(edges, u, v, g, m):
    edges[(u, v)] = g % m
    edges[(v, u)] = inv(g, m)

def undirected_pairs(edges):
    return sorted((u, v) for (u, v) in edges if u < v)

def two_extension(edges, f1, f2, m, new_vertex):
    """Delete directed representatives f1/f2 and insert v -> four endpoints.
    Labels use a=c=0, b=psi(f1), d=psi(f2), hence a^-1 b=f1 and c^-1 d=f2.
    Reverse arcs are inserted automatically.
    """
    out = dict(edges)
    (v1, v2), (v3, v4) = f1, f2
    g1, g2 = out[f1], out[f2]
    assert out[(v2, v1)] == inv(g1, m) and out[(v4, v3)] == inv(g2, m)
    del out[(v1, v2)]; del out[(v2, v1)]
    del out[(v3, v4)]; del out[(v4, v3)]
    add_edge(out, new_vertex, v1, 0, m)
    add_edge(out, new_vertex, v2, g1, m)
    add_edge(out, new_vertex, v3, 0, m)
    add_edge(out, new_vertex, v4, g2, m)
    return out, (g1, g2)

def components(edges, n):
    adj=[[] for _ in range(n)]
    for u,v in undirected_pairs(edges): adj[u].append(v); adj[v].append(u)
    seen=set(); cs=[]
    for s in range(n):
        if s in seen: continue
        q=[s]; seen.add(s); c=[]
        while q:
            u=q.pop(); c.append(u)
            for v in adj[u]:
                if v not in seen: seen.add(v); q.append(v)
        cs.append(c)
    return cs

def diameter(edges, n):
    adj=[[] for _ in range(n)]
    for u,v in undirected_pairs(edges): adj[u].append(v); adj[v].append(u)
    best=0
    for s in range(n):
        d={s:0}; q=deque([s])
        while q:
            u=q.popleft()
            for v in adj[u]:
                if v not in d: d[v]=d[u]+1; q.append(v)
        if len(d)!=n: return float('inf')
        best=max(best,max(d.values()))
    return best

def build_rigid(n=32, m=8):
    # A connected seed, then repeated 2-extensions. Edge count is n+2 here.
    e={}
    for i in range(4): add_edge(e, i, (i+1)%4, (i*3+1)%m, m)
    next_v=4
    while next_v<n:
        pairs=undirected_pairs(e)
        # deterministic, non-identical endpoint pairs when possible
        f1=pairs[(next_v*3) % len(pairs)]
        f2=pairs[(next_v*7+1) % len(pairs)]
        if f1==f2: f2=pairs[(pairs.index(f2)+1)%len(pairs)]
        e,_=two_extension(e, f1, f2, m, next_v); next_v+=1
    return e

def random_graph(n, edge_count, rng):
    allp=[(u,v) for u in range(n) for v in range(u+1,n)]
    chosen=rng.choice(len(allp), edge_count, replace=False)
    e={}
    for j in chosen:
        u,v=allp[j]; add_edge(e,u,v,int(rng.integers(8)),8)
    return e

def reach_fraction(edges, n, steps):
    adj=[set() for _ in range(n)]
    for u,v in undirected_pairs(edges): adj[u].add(v); adj[v].add(u)
    reachable={0}
    for _ in range(steps):
        reachable |= {v for u in reachable for v in adj[u]}
    return len(reachable)/n

def propagation(edges, n, steps=8, d=2):
    # Gain-conditioned linear message passing, row-normalized over undirected neighbors.
    A=np.zeros((n*d,n*d))
    deg=[0]*n
    for u,v in undirected_pairs(edges): deg[u]+=1; deg[v]+=1
    for (u,v),g in edges.items():
        if deg[u]: A[u*d:(u+1)*d,v*d:(v+1)*d]=rotation(g,8)/deg[u]
    x=np.zeros(n*d); x[0]=1.; x[1]=0.
    for _ in range(steps): x=A@x
    return float(np.linalg.norm(x.reshape(n,d),axis=1).sum()), float(np.count_nonzero(np.linalg.norm(x.reshape(n,d),axis=1)>1e-10)) / n

def main():
    n=32; rigid=build_rigid(n); E=len(undirected_pairs(rigid)); rng=np.random.default_rng(SEED)
    random_stats=[]
    for _ in range(100):
        g=random_graph(n,E,rng); random_stats.append((len(components(g,n)),diameter(g,n),reach_fraction(g,n,8),propagation(g,n)[1]))
    de={}
    for u in range(n):
        for v in range(u+1,n): add_edge(de,u,v,(u*13+v*7)%8,8)
    dense={"edges":len(undirected_pairs(de)),"components":len(components(de,n)),"diameter":diameter(de,n),"reach_8_steps":reach_fraction(de,n,8),"active_fraction":propagation(de,n)[1]}
    # Verify every extension's defining equations and reverse-edge inverse law independently.
    check={"reverse_inverse":True,"extension_equations":True}
    e={}; add_edge(e,0,1,3,8); add_edge(e,1,2,5,8); add_edge(e,2,3,1,8); add_edge(e,3,0,6,8)
    for (u,v),g in e.items(): check["reverse_inverse"] &= (e[(v,u)]==inv(g,8))
    old={(u,v):g for (u,v),g in e.items()}
    # Test the general construction with nonzero a,c: b=a*f1 and d=c*f2.
    out=dict(e); del out[(0,1)]; del out[(1,0)]; del out[(2,3)]; del out[(3,2)]
    a,c=3,6
    add_edge(out,4,0,a,8); add_edge(out,4,1,compose(a,old[(0,1)],8),8)
    add_edge(out,4,2,c,8); add_edge(out,4,3,compose(c,old[(2,3)],8),8)
    check["extension_equations"] = (compose(inv(a,8),out[(4,1)],8)==old[(0,1)] and compose(inv(c,8),out[(4,3)],8)==old[(2,3)])
    check["general_extension_reverse_edges"] = all(out[(v,u)]==inv(g,8) for (u,v),g in out.items())
    rs=[len(components(rigid,n)),diameter(rigid,n),reach_fraction(rigid,n,8),propagation(rigid,n)[1]]
    result={"seed":SEED,"n":n,"undirected_edges":E,"math_check":check,
      "gain_rotation_orthogonality_error":float(np.linalg.norm(rotation(3,8).T@rotation(3,8)-np.eye(2))),
      "rigid":{"components":rs[0],"diameter":rs[1],"reach_8_steps":rs[2],"active_fraction":rs[3]},
      "dense_reference":dense,"random_100":{"mean_components":float(np.mean([x[0] for x in random_stats])),"disconnected_rate":float(np.mean([x[0]>1 for x in random_stats])),"mean_diameter_connected":float(np.mean([x[1] for x in random_stats if np.isfinite(x[1])] or [float('inf')])),"mean_reach_8_steps":float(np.mean([x[2] for x in random_stats])),"mean_active_fraction":float(np.mean([x[3] for x in random_stats]))}}
    with open('results.json','w') as f: json.dump(result,f,indent=2)
    print(json.dumps(result,indent=2))

if __name__=='__main__': main()