verify.py
19.1 kB · python · 466 lines
1import cmath2import itertools3import math4import time56# RASTER78ORDERS = 139TWO_PI = 2.0 * math.pi10KIND = {0: "corner", 1: "edge", 2: "corner", 3: "edge", 4: "centre", 5: "edge", 6: "corner", 7: "edge", 8: "corner"}1112def maps(side):13 out = []14 for flip in range(2):15 for turn in range(4):16 m = []17 for f in range(side * side):18 r, c = divmod(f, side)19 if flip:20 r, c = c, r21 for _ in range(turn):22 r, c = c, side - 1 - r23 m.append(r * side + c)24 if m not in out:25 out.append(m)26 return out2728def pair_table(side, group):29 cells = side * side30 index = [-1] * (cells * cells)31 reps = []32 count = 033 for j in range(cells):34 for l in range(j, cells):35 if index[j * cells + l] != -1:36 continue37 for m in group:38 index[m[j] * cells + m[l]] = count39 index[m[l] * cells + m[j]] = count40 reps.append((j, l))41 count += 142 return index, reps, count4344def census(points, index, cells, count):45 out = [0] * count46 for a in range(len(points)):47 row = points[a] * cells48 for b in range(a, len(points)):49 out[index[row + points[b]]] += 150 return out5152def cells_of(bits, cells):53 return [j for j in range(cells) if bits >> j & 1]5455def kronecker(bits):56 base = cells_of(bits, 9)57 return sorted(((p // 3) * 3 + s // 3) * 9 + ((p % 3) * 3 + s % 3) for p in base for s in base)5859def canonical(bits, group):60 best = None61 for m in group:62 x = 063 for j in range(9):64 if bits >> j & 1:65 x |= 1 << m[j]66 best = x if best is None or x < best else best67 return best6869def cube_maps(side):70 axes = [(0, 1, 2), (0, 2, 1), (1, 0, 2), (1, 2, 0), (2, 0, 1), (2, 1, 0)]71 out = []72 for axis in axes:73 for signs in range(8):74 m = []75 for f in range(side ** 3):76 p = (f // (side * side), (f // side) % side, f % side)77 q = [p[axis[0]], p[axis[1]], p[axis[2]]]78 for k in range(3):79 if signs >> k & 1:80 q[k] = side - 1 - q[k]81 m.append((q[0] * side + q[1]) * side + q[2])82 if m not in out:83 out.append(m)84 return out8586def class_count(cells, group):87 index = [-1] * (cells * cells)88 count = 089 for j in range(cells):90 for l in range(j, cells):91 if index[j * cells + l] != -1:92 continue93 for m in group:94 index[m[j] * cells + m[l]] = count95 index[m[l] * cells + m[j]] = count96 count += 197 return count9899def burnside(cells, group):100 total = 0101 for m in group:102 seen = [False] * cells103 cycles = 0104 for j in range(cells):105 if seen[j]:106 continue107 cycles += 1108 x = j109 while not seen[x]:110 seen[x] = True111 x = m[x]112 total += 1 << cycles113 return total // len(group)114115# HARMONICS116117def boxes(side):118 out = []119 for f in range(side * side):120 r, c = divmod(f, side)121 out.append((c / side - 0.5, (c + 1) / side - 0.5, r / side - 0.5, (r + 1) / side - 0.5))122 return out123124def support(shape):125 a1, a2, b1, b2 = shape126 near_x = 0.0 if a1 <= 0.0 <= a2 else min(abs(a1), abs(a2))127 near_y = 0.0 if b1 <= 0.0 <= b2 else min(abs(b1), abs(b2))128 far = max(math.hypot(x, y) for x in (a1, a2) for y in (b1, b2))129 return math.hypot(near_x, near_y), far130131def breaks(side):132 edge = {abs(c / side - 0.5) for c in range(side + 1)} | {0.0}133 reach = math.sqrt(2) / 2134 return sorted({round(math.hypot(a, b), 12) for a in edge for b in edge if math.hypot(a, b) <= reach + 1e-12})135136def arcs(shape, r):137 a1, a2, b1, b2 = shape138 cuts = []139 for a in (a1, a2):140 if abs(a) < r:141 t = math.acos(max(-1.0, min(1.0, a / r)))142 cuts += [t, TWO_PI - t]143 for b in (b1, b2):144 if abs(b) < r:145 s = math.asin(max(-1.0, min(1.0, b / r)))146 cuts += [s % TWO_PI, (math.pi - s) % TWO_PI]147 if not cuts:148 return [(0.0, TWO_PI)] if a1 <= r <= a2 and b1 <= 0.0 <= b2 else []149 cuts.sort()150 out = []151 for i, lo in enumerate(cuts):152 hi = cuts[i + 1] if i + 1 < len(cuts) else cuts[0] + TWO_PI153 if hi - lo < 1e-15:154 continue155 mid = 0.5 * (lo + hi)156 if a1 <= r * math.cos(mid) <= a2 and b1 <= r * math.sin(mid) <= b2:157 out.append((lo, hi))158 return out159160def coefficients(shape, r):161 g = [0j] * ORDERS162 for lo, hi in arcs(shape, r):163 g[0] += hi - lo164 for m in range(1, ORDERS):165 g[m] += (cmath.exp(-1j * m * hi) - cmath.exp(-1j * m * lo)) / (-1j * m)166 return [v / TWO_PI for v in g]167168def quadrature(step, reach):169 nodes = []170 k = 0171 while k * step <= reach:172 for u in ((0.0,) if k == 0 else (k * step, -k * step)):173 s = math.sinh(u) * math.pi / 2174 nodes.append((math.tanh(s), (math.pi / 2) * math.cosh(u) / math.cosh(s) ** 2 * step))175 k += 1176 return nodes177178def gram(side, step, reach, pairs):179 shapes = boxes(side)180 reach_of = [support(shape) for shape in shapes]181 cut = breaks(side)182 nodes = quadrature(step, reach)183 out = [[0.0] * ORDERS for _ in pairs]184 for s in range(len(cut) - 1):185 lo, hi = cut[s], cut[s + 1]186 half, mid = 0.5 * (hi - lo), 0.5 * (hi + lo)187 for x, w in nodes:188 r = mid + half * x189 weight = half * w * TWO_PI * r190 live = {j for j in range(side * side) if reach_of[j][0] <= r <= reach_of[j][1]}191 table = {j: coefficients(shapes[j], r) for j in live}192 for i, (j, l) in enumerate(pairs):193 if j not in live or l not in live:194 continue195 gj, gl = table[j], table[l]196 row = out[i]197 for m in range(ORDERS):198 row[m] += weight * (gj[m] * gl[m].conjugate()).real199 return out, len(nodes), len(cut) - 1200201def spectrum(counts, coefs, mult):202 P = [0.0] * ORDERS203 for c, n in enumerate(counts):204 if n:205 w = mult[c] * n206 row = coefs[c]207 for m in range(ORDERS):208 P[m] += w * row[m]209 return P210211def rank(rows, tolerance):212 work = [row[:] for row in rows]213 width = len(work[0])214 got = 0215 pivots = []216 for col in range(width):217 best = max(range(got, len(work)), key=lambda i: abs(work[i][col]))218 if abs(work[best][col]) < tolerance:219 continue220 work[got], work[best] = work[best], work[got]221 pivots.append(abs(work[got][col]))222 for i in range(len(work)):223 if i == got:224 continue225 factor = work[i][col] / work[got][col]226 for t in range(col, width):227 work[i][t] -= factor * work[got][t]228 got += 1229 return got, pivots230231# CHECKS232233def say(label, value):234 print(" %-58s %s" % (label, value))235236def main():237 clock = time.time()238239 print("block 1: the raster group and its pair classes")240 group = maps(3)241 assert len(group) == 8, len(group)242 index, reps, classes = pair_table(3, group)243 assert classes == 11, classes244 kinds = ["%s-%s" % (KIND[j], KIND[l]) for j, l in reps]245 assert kinds.count("corner-centre") == 1, kinds246 sizes = [len({(min(m[j], m[l]), max(m[j], m[l])) for m in group}) for j, l in reps]247 assert sizes == [4, 8, 4, 4, 8, 2, 4, 4, 4, 2, 1], sizes248 assert sum(sizes) == 45, sizes249 big = maps(9)250 assert len(big) == 8, len(big)251 index2, reps2, classes2 = pair_table(9, big)252 assert classes2 == 461, classes2253 wide = maps(5)254 assert len(wide) == 8 and class_count(25, wide) == 55255 solid_group = cube_maps(3)256 assert len(solid_group) == 48, len(solid_group)257 assert class_count(27, solid_group) == 24258 huge = maps(25)259 wider = class_count(625, huge)260 assert len(huge) == 8 and wider == 24805, wider261 say("pair classes at level 1 and level 2", "%d and %d" % (classes, classes2))262 say("class sizes over the 45 unordered cell pairs", str(sizes))263 say("pair classes on the base-5 plane and the base-3 cube", "55 under 8 symmetries, 24 under 48")264 say("pair classes on the 25 x 25 base-5 level-2 raster", wider)265266 print("block 2: the orbit count from the cycle index and from canonical forms")267 total = burnside(9, group)268 assert total == 102, total269 orbits = {}270 for bits in range(1, 512):271 orbits.setdefault(canonical(bits, group), []).append(bits)272 assert len(orbits) == total - 1 == 101, len(orbits)273 assert sum(len(v) for v in orbits.values()) == 511274 assert max(len(v) for v in orbits.values()) == 8275 say("Burnside average over the eight symmetries", "%d, so %d nonempty" % (total, total - 1))276 say("orbits by canonical form over all 511 codes", "%d, largest %d" % (len(orbits), max(len(v) for v in orbits.values())))277278 print("block 3: the level-1 census takes 97 values on the 101 orbits")279 buckets = {}280 for code in orbits:281 buckets.setdefault(tuple(census(cells_of(code, 9), index, 9, classes)), []).append(code)282 assert len(buckets) == 97, len(buckets)283 tied = sorted(tuple(sorted(v)) for v in buckets.values() if len(v) > 1)284 assert tied == [(45, 105), (61, 121), (78, 102), (94, 118)], tied285 for a, b in tied:286 assert bin(a).count("1") == bin(b).count("1")287 assert canonical(a, group) != canonical(b, group)288 want = {(45, 105): [2, 2, 1, 0, 2, 0, 2, 0, 0, 1, 0], (61, 121): [2, 2, 1, 2, 2, 0, 2, 0, 2, 1, 1], (78, 102): [2, 2, 0, 0, 2, 1, 2, 1, 0, 0, 0], (94, 118): [2, 2, 0, 2, 2, 1, 2, 1, 2, 0, 1]}289 for (a, b), counts in want.items():290 assert census(cells_of(a, 9), index, 9, classes) == census(cells_of(b, 9), index, 9, classes) == counts, (a, b)291 assert (61, 121) == (45 | 16, 105 | 16) and (94, 118) == (78 | 16, 102 | 16)292 drawn = {45: [(0, 0), (2, 0), (0, 1), (2, 1)], 105: [(0, 0), (0, 1), (2, 1), (0, 2)]}293 for code, filled in drawn.items():294 assert sum(1 << (r * 3 + c) for c, r in filled) == code, code295 say("distinct level-1 pair censuses on the 101 orbits", len(buckets))296 say("homometric pairs, equal in all 11 class counts", str(tied))297 say("the census shared by 45 and 105", str(want[(45, 105)]))298 say("the second two pairs are the first two plus the centre cell", "61 = 45 + centre, 94 = 78 + centre")299 say("the column and row of each cell drawn in figure 1", "45 at %s, 105 at %s" % (drawn[45], drawn[105]))300301 print("block 4: the level-2 census separates all four homometric pairs")302 for a, b in tied:303 first = census(kronecker(a), index2, 81, classes2)304 second = census(kronecker(b), index2, 81, classes2)305 moved = sum(1 for x, y in zip(first, second) if x != y)306 assert moved > 0, (a, b)307 say("%d against %d: level-2 class counts that differ" % (a, b), "%d of %d" % (moved, classes2))308309 print("block 5: the level-1 Gram matrix is constant on the pair classes")310 every = [(j, k) for j in range(9) for k in range(9)]311 flat, nodes, segments = gram(3, 0.06, 3.4, every)312 Q = [[[flat[j * 9 + k][m] for k in range(9)] for j in range(9)] for m in range(ORDERS)]313 scale = max(abs(Q[m][j][k]) for m in range(ORDERS) for j in range(9) for k in range(9))314 spread = 0.0315 for m in range(ORDERS):316 seen = [[] for _ in range(classes)]317 for j in range(9):318 for k in range(9):319 seen[index[j * 9 + k]].append(Q[m][j][k])320 spread = max(spread, max(max(v) - min(v) for v in seen))321 assert spread < 1e-15, spread322 coarse, _, _ = gram(3, 0.03, 4.0, every)323 drift = max(abs(flat[i][m] - coarse[i][m]) for i in range(81) for m in range(ORDERS))324 assert drift < 1e-15, drift325 say("radial segments and quadrature nodes on each", "%d and %d" % (segments, nodes))326 say("worst spread of Q_m inside one pair class", "%.2e on scale %.3e" % (spread, scale))327 say("worst drift against a finer quadrature rule", "%.2e" % drift)328329 print("block 6: the two exact relations that cost the census two directions")330 hole = max(abs(Q[m][0][4]) for m in range(ORDERS))331 assert hole == 0.0, hole332 similar = max(abs(sum(Q[m][j][k] for j in range(9) for k in range(9)) - 9 * Q[m][4][4]) for m in range(ORDERS))333 assert similar < 1e-15, similar334 ordered = [0] * classes335 for j in range(9):336 for k in range(9):337 ordered[index[j * 9 + k]] += 1338 assert ordered == [4, 16, 8, 8, 16, 4, 4, 8, 8, 4, 1], ordered339 assert sum(ordered) == 81340 byclass = max(abs(sum(ordered[c] * Q[m][reps[c][0]][reps[c][1]] for c in range(classes)) - 9 * Q[m][4][4]) for m in range(ORDERS))341 assert byclass < 1e-15, byclass342 say("Q_m at a corner against the centre, every order", "%.1e" % hole)343 say("ordered pairs in each class", str(ordered))344 say("worst |P_m(square) - 9 P_m(centre cell)|", "%.2e cell by cell, %.2e by class" % (similar, byclass))345346 print("block 7: the thirteen orders reach rank 9, even 6, odd 3")347 C = [[Q[m][reps[c][0]][reps[c][1]] for c in range(classes)] for m in range(ORDERS)]348 for tolerance in (1e-9, 1e-11, 1e-13):349 got, pivots = rank(C, tolerance)350 assert got == 9, (tolerance, got)351 even, _ = rank([C[m] for m in range(0, ORDERS, 2)], 1e-11)352 odd, _ = rank([C[m] for m in range(1, ORDERS, 2)], 1e-11)353 assert (even, odd) == (6, 3), (even, odd)354 say("rank of the 13 orders against the 11 classes", "%d, stable from 1e-9 to 1e-13" % got)355 say("smallest pivot against the largest coefficient", "%.3e against %.3e" % (min(pivots), scale))356 say("even orders and odd orders", "%d and %d" % (even, odd))357358 print("block 8: the half-turn involution caps the odd orders at 3")359 half = [8 - j for j in range(9)]360 assert half in group, half361 tau = [index[half[reps[c][0]] * 9 + reps[c][1]] for c in range(classes)]362 assert all(tau[tau[c]] == c for c in range(classes)), tau363 fixed = sum(1 for c in range(classes) if tau[c] == c)364 assert fixed == 5, fixed365 assert (classes - fixed) // 2 == 3366 for m in range(ORDERS):367 sign = 1 if m % 2 == 0 else -1368 worst = max(abs(C[m][tau[c]] - sign * C[m][c]) for c in range(classes))369 assert worst < 1e-15, (m, worst)370 say("classes fixed by half-turning one member", "%d of %d" % (fixed, classes))371 say("antisymmetric dimension, the cap on the odd orders", (classes - fixed) // 2)372373 print("block 9: the census route and the cell route agree, and the homometric pairs share every P_m")374 mult = [1 if reps[c][0] == reps[c][1] else 2 for c in range(classes)]375 coefs = [[Q[m][reps[c][0]][reps[c][1]] for m in range(ORDERS)] for c in range(classes)]376 level1 = {}377 routes = 0.0378 for bits in range(1, 512):379 points = cells_of(bits, 9)380 direct = [sum(Q[m][j][k] for j in points for k in points) for m in range(ORDERS)]381 P = spectrum(census(points, index, 9, classes), coefs, mult)382 routes = max(routes, max(abs(x - y) for x, y in zip(direct, P)))383 level1[bits] = direct384 assert routes < 1e-15, routes385 say("worst gap between the two routes over all 511 codes", "%.2e" % routes)386 for a, b in tied:387 gap = max(abs(x - y) for x, y in zip(level1[a], level1[b]))388 size = max(abs(x) for x in level1[a])389 assert gap < 1e-15, (a, b, gap)390 say("%d against %d: worst |P_m| gap" % (a, b), "%.2e on scale %.3e" % (gap, size))391 spectra = {}392 for bits in range(1, 512):393 spectra.setdefault(tuple(round(x, 9) for x in level1[bits]), []).append(bits)394 assert len(spectra) == 97, len(spectra)395 say("distinct level-1 spectra over all 511 codes at 1e-9", len(spectra))396397 print("block 10: the level-2 spectrum splits the 511 codes into 101")398 mult2 = [1 if reps2[c][0] == reps2[c][1] else 2 for c in range(classes2)]399 coefs2, nodes2, segments2 = gram(9, 0.06, 3.4, reps2)400 gens = (big[1], big[4])401 seen = {tuple(range(81))}402 edge = list(seen)403 while edge:404 cur = edge.pop()405 for g in gens:406 nxt = tuple(g[j] for j in cur)407 if nxt not in seen:408 seen.add(nxt)409 edge.append(nxt)410 assert len(seen) == 8 and all(list(g) in big for g in seen), len(seen)411 turned = 0.0412 for g in gens:413 check2, _, _ = gram(9, 0.06, 3.4, [(g[j], g[l]) for j, l in reps2])414 turned = max(turned, max(abs(coefs2[c][m] - check2[c][m]) for c in range(classes2) for m in range(ORDERS)))415 assert turned < 1e-15, turned416 scale2 = max(abs(coefs2[c][m]) for c in range(classes2) for m in range(ORDERS))417 finer2, _, _ = gram(9, 0.03, 4.0, reps2)418 drift2 = max(abs(coefs2[c][m] - finer2[c][m]) for c in range(classes2) for m in range(ORDERS))419 assert drift2 < 1e-15, drift2420 level2 = {}421 for bits in range(1, 512):422 level2[bits] = spectrum(census(kronecker(bits), index2, 81, classes2), coefs2, mult2)423 grouped = {}424 for bits in range(1, 512):425 grouped.setdefault(tuple(round(x, 9) for x in level2[bits]), []).append(bits)426 assert len(grouped) == 101, len(grouped)427 assert max(len(v) for v in grouped.values()) == 8428 for v in grouped.values():429 assert len({canonical(b, group) for b in v}) == 1, v430 closest = None431 for a, b in itertools.combinations(sorted(orbits), 2):432 gap = max(abs(x - y) for x, y in zip(level2[a], level2[b]))433 rel = gap / max(max(map(abs, level2[a])), max(map(abs, level2[b])))434 if closest is None or rel < closest[0]:435 closest = (rel, gap, a, b)436 assert closest[0] > 1e-3, closest437 say("radial segments and quadrature nodes at level 2", "%d and %d" % (segments2, nodes2))438 say("worst drift against a finer quadrature rule at level 2", "%.3e on scale %.3e" % (drift2, scale2))439 say("worst drift of Q_m under the two generators of the group", "%.2e" % turned)440 say("distinct level-2 spectra over all 511 codes", "%d, largest bucket %d" % (len(grouped), max(len(v) for v in grouped.values())))441 say("every bucket is one symmetry orbit", "yes, %d orbits" % len(orbits))442 say("closest two orbits, codes %d and %d" % (closest[2], closest[3]), "relative %.3e, absolute %.3e" % (closest[0], closest[1]))443444 print("block 11: the truncation is a truncation")445 share = {}446 for bits in range(1, 512):447 energy = bin(bits).count("1") / 9.0448 share[bits] = (level1[bits][0] + 2 * sum(level1[bits][m] for m in range(1, ORDERS))) / energy449 solid = 511450 caught = share[solid]451 assert 0.977 < caught < 0.978, caught452 best = max(share.values())453 worst = min(share.values())454 top = sorted(b for b in share if share[b] > best - 1e-12)455 low = sorted(b for b in share if share[b] < worst + 1e-12)456 assert top == [16, 511], top457 assert low == [1, 4, 64, 256], low458 assert abs(best - caught) < 1e-12, (best, caught)459 say("share of the solid square's angular energy at m <= 12", "%.6f of 1" % caught)460 say("the largest share over all 511 codes and where", "%.6f at codes %s" % (best, top))461 say("the smallest share over all 511 codes and where", "%.6f at codes %s" % (worst, low))462463 print("all green in %.1f seconds" % (time.time() - clock))464465if __name__ == "__main__":466 main()