verify.py
12.6 kB · python · 352 lines
1import math2from fractions import Fraction3from itertools import permutations, product45# COPRIME67DESIGNS = [8 ("or-triangle", 2, 2, [(1, 0), (0, 1), (1, 1)], 11,9 [3, 6, 22, 58, 200, 576, 1798, 5174, 15944, 47744, 143808]),10 ("gasket", 2, 2, [(0, 0), (1, 0), (0, 1)], 11,11 [2, 4, 12, 34, 122, 362, 1130, 3406, 10506, 31550, 95260]),12 ("carpet", 3, 2, [(0, 0), (0, 1), (0, 2), (1, 0), (1, 2), (2, 0), (2, 1), (2, 2)], 5,13 [4, 32, 274, 2320, 19178]),14 ("vicsek-plus", 3, 2, [(1, 0), (0, 1), (1, 1), (2, 1), (1, 2)], 7,15 [5, 16, 90, 418, 2178, 10560, 54120]),16 ("q4-rows", 4, 2, [(0, 0), (0, 1), (0, 2), (0, 3), (1, 0), (1, 1), (1, 2), (1, 3)], 5,17 [5, 37, 302, 2340, 19100]),18 ("q4-border", 4, 2, [(0, 0), (0, 1), (0, 3), (1, 0), (1, 3), (3, 0), (3, 2), (3, 3)], 5,19 [4, 33, 318, 2690, 22675]),20 ("q5-binary-3d", 5, 3, [(0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0), (1, 0, 1), (1, 1, 0), (1, 1, 1)], 6,21 [7, 42, 294, 2088, 14412, 101346]),22 ("q5-axes-3d", 5, 3, [(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1), (2, 0, 0), (0, 2, 0), (0, 0, 2)], 6,23 [3, 21, 189, 1467, 10911, 79425]),24 ("q6-first", 6, 2, [(0, 0), (0, 1), (1, 0), (1, 1), (2, 3), (3, 2), (4, 5), (5, 4)], 5,25 [7, 43, 403, 3105, 25881]),26 ("q6-second", 6, 2, [(0, 0), (0, 3), (1, 1), (1, 2), (2, 0), (2, 1), (4, 0), (5, 5)], 6,27 [3, 25, 230, 1789, 14918, 118380]),28]2930SPONGE = [v for v in product(range(3), repeat=3) if sum(1 for c in v if c == 1) <= 1]3132def die(what, got, want):33 raise AssertionError("%s: got %r, want %r" % (what, got, want))3435def check(what, got, want):36 if got != want:37 die(what, got, want)3839def close(what, got, want, tol):40 if abs(got - want) > tol:41 die("%s (tol %g)" % (what, tol), got, want)4243def zeta3():44 n = 10000045 total = sum(1.0 / j ** 3 for j in range(1, n + 1))46 return total + 1.0 / (2 * n ** 2) - 1.0 / (2 * n ** 3) + 1.0 / (4 * n ** 4)4748def primes_of(m):49 out = []50 d = 251 while d * d <= m:52 if m % d == 0:53 out.append(d)54 while m % d == 0:55 m //= d56 d += 157 if m > 1:58 out.append(m)59 return out6061def bracket(q, F):62 k = len(F)63 ps = primes_of(q)64 total = Fraction(0)65 for mask in range(1 << len(ps)):66 e = 167 sign = 168 for i, p in enumerate(ps):69 if mask >> i & 1:70 e *= p71 sign = -sign72 ke = sum(1 for v in F if all(c % e == 0 for c in v))73 total += sign * Fraction(ke, k)74 return total7576def predicted(q, D, F, z):77 val = float(bracket(q, F)) / z78 for p in primes_of(q):79 val /= 1.0 - p ** (-D)80 return val8182def gcd_of(x):83 g = 084 for c in x:85 g = math.gcd(g, c)86 return g8788def level(q, F, D, n):89 pts = [tuple([0] * D)]90 for _ in range(n):91 pts = [tuple(q * a + b for a, b in zip(x, f)) for x in pts for f in F]92 return pts9394def walk(q, F, D, n):95 pts = [tuple([0] * D)]96 for _ in range(n):97 pts = [tuple(q * a + b for a, b in zip(x, f)) for x in pts for f in F]98 yield pts99100# CONSTANTS101102def constants(z2, z3):103 named = [104 ("gasket", 2, 2, [(0, 0), (1, 0), (0, 1)], Fraction(2, 3), 16.0 / (3 * math.pi ** 2)),105 ("or-triangle", 2, 2, [(1, 0), (0, 1), (1, 1)], Fraction(1), 8.0 / math.pi ** 2),106 ("carpet", 3, 2, [(0, 0), (0, 1), (0, 2), (1, 0), (1, 2), (2, 0), (2, 1), (2, 2)],107 Fraction(7, 8), 189.0 / (32 * math.pi ** 2)),108 ("vicsek-plus", 3, 2, [(1, 0), (0, 1), (1, 1), (2, 1), (1, 2)],109 Fraction(1), 27.0 / (4 * math.pi ** 2)),110 ("q6-second", 6, 2, [(0, 0), (0, 3), (1, 1), (1, 2), (2, 0), (2, 1), (4, 0), (5, 5)],111 Fraction(1, 2), 9.0 / (2 * math.pi ** 2)),112 ("q4-rows", 4, 2, DESIGNS[4][3], Fraction(3, 4), 6.0 / math.pi ** 2),113 ("q4-border", 4, 2, DESIGNS[5][3], Fraction(7, 8), 7.0 / math.pi ** 2),114 ("q5-binary-3d", 5, 3, DESIGNS[6][3], Fraction(1), (125.0 / 124.0) / z3),115 ("q5-axes-3d", 5, 3, DESIGNS[7][3], Fraction(6, 7), (375.0 / 434.0) / z3),116 ("q6-first", 6, 2, DESIGNS[8][3], Fraction(7, 8), 63.0 / (8 * math.pi ** 2)),117 ]118 for name, q, D, F, wantB, wantD in named:119 check("bracket %s" % name, bracket(q, F), wantB)120 close("delta %s" % name, predicted(q, D, F, z2 if D == 2 else z3), wantD, 1e-12)121 print("constant %s: B = %s, delta = %.12f" % (name, wantB, wantD))122 check("bracket sponge", bracket(3, SPONGE), Fraction(19, 20))123 sponge = (513.0 / 520.0) / z3124 close("delta sponge", predicted(3, 3, SPONGE, z3), sponge, 1e-12)125 print("constant sponge: B = 19/20, delta = %.12f" % sponge)126 print("brackets at composite bases separate equal k: 3/4 vs 7/8 at q = 4, 7/8 vs 1/2 at q = 6")127128# ENUMERATION129130def enumeration(z2, z3):131 rows = 0132 points = 0133 for name, q, D, F, N, want in DESIGNS:134 k = len(F)135 B = bracket(q, F)136 check("%s lattice index" % name, lattice_index(F), 1)137 got = []138 for n, pts in enumerate(walk(q, F, D, N), start=1):139 points += len(pts)140 a = 0141 coprime_to_q = 0142 for x in pts:143 g = gcd_of(x)144 if g == 1:145 a += 1146 if math.gcd(g, q) == 1:147 coprime_to_q += 1148 check("%s A(%d)" % (name, n), a, want[n - 1])149 check("%s base identity at n = %d" % (name, n), Fraction(coprime_to_q), B * k ** n)150 got.append(a)151 rows += 1152 d = predicted(q, D, F, z2 if D == 2 else z3)153 print("%-13s q=%d D=%d k=%2d B=%-4s A(%d)=%d err %+.6f"154 % (name, q, D, k, B, N, got[-1], got[-1] / float(k ** N) - d))155 check("level rows", rows, 67)156 print("enumerated %d level rows, %d points, no sampling" % (rows, points))157 print("all ten designs have lattice index 1, so all ten are spanning")158159# FACTOR160161def factor():162 q, D = 6, 2163 F = DESIGNS[9][3]164 k = len(F)165 check("k_2", sum(1 for v in F if all(c % 2 == 0 for c in v)), 3)166 check("k_3", sum(1 for v in F if all(c % 3 == 0 for c in v)), 2)167 check("k_6", sum(1 for v in F if all(c % 6 == 0 for c in v)), 1)168 check("bracket", bracket(q, F), Fraction(1, 2))169 for n, pts in enumerate(walk(q, F, D, 6), start=1):170 got = sum(1 for x in pts if math.gcd(gcd_of(x), 6) == 1)171 check("gcd coprime to 6 at n = %d" % n, got, k ** n // 2)172 if n >= 3:173 naive = 0.46875 * k ** n174 if abs(got - naive) < 1:175 die("naive product at n = %d" % n, naive, "a value differing from %d" % got)176 print("n=%d |S_n|=%-7d coprime to 6 = %-7d naive = %.0f" % (n, k ** n, got, naive))177 print("the bracket 1/2 is not (1-3/8)(1-2/8) = 0.46875")178179# SHARPNESS180181def sharpness(z2):182 dust = [(0, 0), (0, 2), (2, 0), (2, 2)]183 for n, pts in enumerate(walk(3, dust, 2, 8), start=1):184 check("cantor dust A(%d)" % n, sum(1 for x in pts if gcd_of(x) == 1), 0)185 check("dust bracket", bracket(3, dust), Fraction(3, 4))186 close("dust predicted delta", predicted(3, 2, dust, z2), 81.0 / (16 * math.pi ** 2), 1e-12)187 check("dust lattice index", lattice_index(dust), 4)188 print("cantor dust {0,2}^2 at base 3: A(n) = 0 for 1 <= n <= 8, k = 4 > 3, predicted 0.512938")189 sheared = [(0, 0), (0, 1), (1, 0), (1, 1)]190 check("sheared bracket", bracket(3, sheared), Fraction(3, 4))191 want = 81.0 / (16 * math.pi ** 2)192 close("sheared delta", predicted(3, 2, sheared, z2), want, 1e-12)193 for n, pts in enumerate(walk(3, dust, 2, 8), start=1):194 got = sum(1 for x in pts if gcd_of(x) == 2)195 base = sum(1 for x in level(3, sheared, 2, n) if gcd_of(x) == 1)196 check("shear bijection at n = %d" % n, got, base)197 print("sheared design {0,1}^2: gcd = 2 density -> 81/(16 pi^2) = %.12f" % want)198199# BOX200201def box():202 trials = [("gasket", 2, 2, DESIGNS[1][3]), ("carpet", 3, 2, DESIGNS[2][3]),203 ("q6-second", 6, 2, DESIGNS[9][3])]204 tested = 0205 for name, q, D, F in trials:206 k = len(F)207 alpha = math.log(k) / math.log(q)208 for n, pts in enumerate(walk(q, F, D, 6), start=1):209 tally = {}210 for x in pts:211 g = gcd_of(x)212 if g:213 tally[g] = tally.get(g, 0) + 1214 for m in range(2, 61):215 got = sum(c for g, c in tally.items() if g % m == 0)216 bound = (q + 1) ** D * k ** n * m ** (-alpha)217 if got > bound:218 die("box bound %s n=%d m=%d" % (name, n, m), got, "at most %.4f" % bound)219 tested += 1220 print("box bound holds in %d exact cases, m <= 60, n <= 6" % tested)221222# CENSUS223224CUBE = list(product([0, 1], repeat=4))225CIDX = {v: i for i, v in enumerate(CUBE)}226227def generators():228 out = []229 for a in range(3):230 p = list(range(4))231 p[a], p[a + 1] = p[a + 1], p[a]232 out.append(tuple(CIDX[tuple(v[p[j]] for j in range(4))] for v in CUBE))233 out.append(tuple(CIDX[(v[0] ^ 1, v[1], v[2], v[3])] for v in CUBE))234 return out235236def burnside(dim):237 verts = list(product([0, 1], repeat=dim))238 idx = {v: i for i, v in enumerate(verts)}239 order = 0240 total = 0241 for perm in permutations(range(dim)):242 for flip in product([0, 1], repeat=dim):243 order += 1244 img = [idx[tuple(v[perm[j]] ^ flip[j] for j in range(dim))] for v in verts]245 seen = [False] * len(verts)246 cycles = 0247 for i in range(len(verts)):248 if not seen[i]:249 cycles += 1250 j = i251 while not seen[j]:252 seen[j] = True253 j = img[j]254 total += 2 ** cycles255 return order, total // order256257def lattice_index(F):258 v0 = F[0]259 rows = [[a - b for a, b in zip(v, v0)] for v in F[1:]]260 dim = len(v0)261 r = 0262 for c in range(dim):263 pivot = None264 for i in range(r, len(rows)):265 if rows[i][c]:266 pivot = i267 break268 if pivot is None:269 return 0270 rows[r], rows[pivot] = rows[pivot], rows[r]271 for i in range(r + 1, len(rows)):272 while rows[i][c]:273 f = rows[r][c] // rows[i][c]274 rows[r] = [a - f * b for a, b in zip(rows[r], rows[i])]275 rows[r], rows[i] = rows[i], rows[r]276 r += 1277 det = 1278 for c in range(dim):279 det *= rows[c][c]280 return abs(det)281282def census():283 for dim, want in [(0, 2), (1, 3), (2, 6), (3, 22)]:284 order, got = burnside(dim)285 check("orbits at D = %d" % dim, got, want)286 parent = list(range(65536))287288 def find(x):289 while parent[x] != x:290 parent[x] = parent[parent[x]]291 x = parent[x]292 return x293294 gens = generators()295 for mask in range(65536):296 a = find(mask)297 for g in gens:298 image = 0299 for i in range(16):300 if mask >> i & 1:301 image |= 1 << g[i]302 b = find(image)303 if a != b:304 lo, hi = (a, b) if a < b else (b, a)305 parent[hi] = lo306 a = lo307 reps = sorted({find(mask) for mask in range(65536)})308 check("orbits at D = 4", len(reps), 402)309 k_ge_2 = k_gt_2 = spanning = 0310 signatures = {}311 for mask in reps:312 F = [CUBE[i] for i in range(16) if mask >> i & 1]313 k = len(F)314 if k >= 2:315 k_ge_2 += 1316 if k > 2:317 k_gt_2 += 1318 if lattice_index(F) == 1:319 spanning += 1320 if k < 2:321 continue322 pts = [(0, 0, 0, 0)]323 sig = []324 for _ in range(4):325 pts = [(2 * x[0] + f[0], 2 * x[1] + f[1], 2 * x[2] + f[2], 2 * x[3] + f[3])326 for x in pts for f in F]327 sig.append(sum(1 for x in pts if gcd_of(x) == 1))328 signatures.setdefault(tuple(sig), []).append(mask)329 check("representatives with k >= 2", k_ge_2, 400)330 check("representatives with k > 2", k_gt_2, 396)331 check("spanning with k > 2", spanning, 336)332 check("orbits outside the sufficient condition", 402 - spanning, 66)333 check("distinct signatures", len(signatures), 189)334 check("collision groups", sum(1 for v in signatures.values() if len(v) > 1), 87)335 print("D = 4 base 2: 402 orbits, 400 with k >= 2, 396 with k > 2, 336 spanning and k > 2")336 print("400 eligible representatives realize 189 distinct (A(1),A(2),A(3),A(4)), 87 repeated")337 print("orbit counts for D = 0..4: 2, 3, 6, 22, 402")338339def main():340 z2 = math.pi ** 2 / 6.0341 z3 = zeta3()342 close("zeta(3)", z3, 1.2020569031595942854, 1e-11)343 constants(z2, z3)344 enumeration(z2, z3)345 factor()346 sharpness(z2)347 box()348 census()349 print("all green")350351if __name__ == "__main__":352 main()