node_stack.py
16.0 kB · python · 449 lines
1import hashlib2from fractions import Fraction3from math import floor, gcd, lcm4from pathlib import Path56import numpy as np7from PIL import Image89HERE = Path(__file__).resolve().parent10PARITY_CHECK_N = 1511PARITY_RENDER_N = 3112CARPET_CHECKS = ((12, 1), (6, 2))13CARPET_CORNER_N = 1214FAREY_MAX_N = 3015DISCREPANCY_QS = (30, 60, 90)16TOP = 617RENDER_R = 102418SAVE_R = 51219SHEET_R = 25620CARPET_RENDER_L = 221CARPET_RENDER_N = 1222INK_FLOOR = 702324def parity_inked(i, j, n):25 return 0 <= i < n and 0 <= j < n and i % 2 == 1 and j % 2 == 12627def parity_edges(n):28 v = set((k, j) for k in range(n + 1) for j in range(n)29 if parity_inked(k - 1, j, n) != parity_inked(k, j, n))30 h = set((i, m) for m in range(n + 1) for i in range(n)31 if parity_inked(i, m - 1, n) != parity_inked(i, m, n))32 return v, h3334def parity_edges_form(n):35 v = set((k, j) for k in range(1, n) for j in range(1, n, 2))36 h = set((i, m) for m in range(1, n) for i in range(1, n, 2))37 return v, h3839def parity_corners(n):40 return set((Fraction(k, n), Fraction(m, n))41 for i in range(1, n, 2) for j in range(1, n, 2)42 for k in (i, i + 1) for m in (j, j + 1))4344def parity_corners_form(n):45 return set((Fraction(k, n), Fraction(m, n))46 for k in range(1, n) for m in range(1, n))4748def parity_corner_stack(n_max):49 out = {}50 for n in range(1, n_max + 1, 2):51 for p in parity_corners(n):52 out[p] = out.get(p, 0) + 153 return out5455def parity_corner_form(b, d, n_max):56 m = lcm(b, d)57 if m % 2 == 0:58 return 059 return (n_max // m + 1) // 26061def parity_corner_predicted(n_max):62 out = {}63 for b in range(1, n_max + 1, 2):64 for d in range(1, n_max + 1, 2):65 if lcm(b, d) > n_max:66 continue67 bright = parity_corner_form(b, d, n_max)68 for a in range(1, b):69 if gcd(a, b) != 1:70 continue71 for c in range(1, d):72 if gcd(c, d) == 1:73 out[(Fraction(a, b), Fraction(c, d))] = bright74 return out7576def carpet_inked(i, j, level):77 for _ in range(level):78 if i % 3 == 1 and j % 3 == 1:79 return False80 i //= 381 j //= 382 return True8384def carpet_cell(k, j, level, side):85 s = 3 ** level86 return 0 <= k < side and 0 <= j < side and carpet_inked(k % s, j % s, level)8788def carpet_edges(n, level):89 side = 3 ** level * n90 v = set((k, j) for k in range(side + 1) for j in range(side)91 if carpet_cell(k - 1, j, level, side) != carpet_cell(k, j, level, side))92 h = set((i, m) for m in range(side + 1) for i in range(side)93 if carpet_cell(i, m - 1, level, side) != carpet_cell(i, m, level, side))94 return v, h9596def edge_residues(level):97 s = 3 ** level98 rows = {}99 for j in range(s):100 r = [x for x in range(s)101 if carpet_inked((j - 1) % s, x, level) != carpet_inked(j, x, level)]102 if r:103 rows[j] = r104 return sorted(rows), rows105106def carpet_edge_lines(n, level):107 side = 3 ** level * n108 v, _ = carpet_edges(n, level)109 return set(Fraction(k, side) for (k, j) in v if 0 < k < side)110111def carpet_line_lit(a, b, n, level, residues):112 s = 3 ** level113 if (s * n) % b:114 return False115 return (a * (s * n // b)) % s in residues116117def v3(x):118 t = 0119 while x % 3 == 0:120 x //= 3121 t += 1122 return t123124def carpet_lines_reach(n_max, level):125 out = set()126 for b in range(3, 3 ** level * n_max + 1, 3):127 if b // 3 ** min(v3(b), level) > n_max:128 continue129 for a in range(1, b):130 if gcd(a, b) == 1:131 out.add(Fraction(a, b))132 return out133134def carpet_edge_segments(n_max, level):135 side_lines = {}136 for n in range(1, n_max + 1):137 side = 3 ** level * n138 v, _ = carpet_edges(n, level)139 for (k, j) in v:140 if k == 0 or k == side:141 continue142 side_lines.setdefault(Fraction(k, side), []).append(143 (Fraction(j, side), Fraction(j + 1, side)))144 return side_lines145146def top_segments(lines, count):147 out = []148 for x, ivs in lines.items():149 pts = sorted(set(p for iv in ivs for p in iv))150 for i in range(len(pts) - 1):151 lo, hi = pts[i], pts[i + 1]152 b = sum(1 for (u, w) in ivs if u <= lo and hi <= w)153 if b:154 out.append((b, hi - lo, x, lo, hi))155 out.sort(key=lambda t: (-t[0], -t[1], t[2], t[3]))156 return out[:count]157158def parity_edge_lines(n):159 v, _ = parity_edges(n)160 return set(Fraction(k, n) for (k, j) in v)161162def parity_line_brightness(b, n_max):163 if b % 2 == 0:164 return 0165 return (n_max // b + 1) // 2166167def carpet_line_brightness(b, n_max, level):168 s = v3(b)169 if s == 0:170 return 0171 return n_max // (b // 3 ** min(s, level)) - n_max // b172173def check_line_brightness_parity(n_max):174 count = {}175 for n in range(1, n_max + 1, 2):176 for f in parity_edge_lines(n):177 count[f] = count.get(f, 0) + 1178 bad = sum(1 for f, c in count.items() if parity_line_brightness(f.denominator, n_max) != c)179 invented = sum(1 for b in range(1, n_max + 1) for a in range(1, b)180 if gcd(a, b) == 1 and parity_line_brightness(b, n_max) != count.get(Fraction(a, b), 0))181 return len(count), bad + invented182183def check_line_brightness_carpet(n_max, level):184 count = {}185 for n in range(1, n_max + 1):186 for f in carpet_edge_lines(n, level):187 count[f] = count.get(f, 0) + 1188 bad = sum(1 for f, c in count.items()189 if carpet_line_brightness(f.denominator, n_max, level) != c)190 invented = sum(1 for b in range(1, 3 ** level * n_max + 1) for a in range(1, b)191 if gcd(a, b) == 1192 and carpet_line_brightness(b, n_max, level) != count.get(Fraction(a, b), 0))193 return len(count), bad + invented194195def parity_segment_check(n_max):196 odds = list(range(1, n_max + 1, 2))197 lines = {}198 for n in odds:199 v, _ = parity_edges(n)200 for (k, j) in v:201 lines.setdefault(Fraction(k, n), []).append((n, Fraction(j, n), Fraction(j + 1, n)))202 pts = sorted(set(Fraction(j, n) for n in odds for j in range(n + 1)))203 mids = [(pts[i] + pts[i + 1]) / 2 for i in range(len(pts) - 1)]204 bad = 0205 for x, ivs in lines.items():206 b = x.denominator207 for y in mids:208 literal = sum(1 for (n, u, w) in ivs if u <= y <= w)209 form = sum(1 for n in odds if n % b == 0 and floor(n * y) % 2 == 1)210 bad += literal != form211 return len(lines) * len(mids), bad212213def carpet_segment_check(n_max, level):214 s = 3 ** level215 _, rows = edge_residues(level)216 lines = {}217 for n in range(1, n_max + 1):218 side = s * n219 v, _ = carpet_edges(n, level)220 for (k, j) in v:221 if 0 < k < side:222 lines.setdefault(Fraction(k, side), []).append(223 (n, Fraction(j, side), Fraction(j + 1, side)))224 pts = sorted(set(Fraction(j, s * n) for n in range(1, n_max + 1) for j in range(s * n + 1)))225 mids = [(pts[i] + pts[i + 1]) / 2 for i in range(len(pts) - 1)]226 bad = 0227 for x, ivs in lines.items():228 a, b = x.numerator, x.denominator229 for y in mids:230 literal = sum(1 for (n, u, w) in ivs if u <= y <= w)231 form = 0232 for n in range(1, n_max + 1):233 side = s * n234 if side % b:235 continue236 j = (a * (side // b)) % s237 if j in rows and floor(side * y) % s in rows[j]:238 form += 1239 bad += literal != form240 return len(lines) * len(mids), bad241242def carpet_corners(n, level):243 side = 3 ** level * n244 out = set()245 for k in range(side + 1):246 for m in range(side + 1):247 if any(carpet_cell(k - 1 + du, m - 1 + dv, level, side)248 for du in (0, 1) for dv in (0, 1)):249 out.add((Fraction(k, side), Fraction(m, side)))250 return out251252def carpet_corner_stack(n_max, level):253 out = {}254 for n in range(1, n_max + 1):255 for p in carpet_corners(n, level):256 out[p] = out.get(p, 0) + 1257 return out258259def carpet_corner_form(b, d, n_max):260 m = lcm(b, d)261 return n_max // (m // gcd(m, 3))262263def totients(q):264 phi = list(range(q + 1))265 for p in range(2, q + 1):266 if phi[p] == p:267 for k in range(p, q + 1, p):268 phi[k] -= phi[k] // p269 return phi270271def farey_sizes(q):272 phi = totients(q)273 return sum(phi[1:]), sum(phi[k] for k in range(3, q + 1, 3))274275def farey_list(q, step=1):276 out = []277 for b in range(step, q + 1, step):278 for a in range(1, b + 1):279 if gcd(a, b) == 1:280 out.append(Fraction(a, b))281 out.sort()282 return out283284def landau(nodes):285 m = len(nodes)286 total = Fraction(0)287 for i, f in enumerate(nodes, start=1):288 total += abs(f - Fraction(i, m))289 return total290291def digest(lines):292 return hashlib.sha256("\n".join(lines).encode("utf-8")).hexdigest()293294def stack_digest(stack):295 return digest(["%s,%s:%d" % (p[0], p[1], b) for p, b in sorted(stack.items())])296297def check_parity_edges(n_max):298 bad = 0299 for n in range(1, n_max + 1, 2):300 if parity_edges(n) != parity_edges_form(n):301 bad += 1302 return bad303304def check_parity_corners(n_max):305 bad = 0306 for n in range(1, n_max + 1, 2):307 if parity_corners(n) != parity_corners_form(n):308 bad += 1309 return bad310311def check_carpet_lines(n_max, level):312 residues, _ = edge_residues(level)313 rs = set(residues)314 bad = 0315 for n in range(1, n_max + 1):316 literal = carpet_edge_lines(n, level)317 form = set(f for f in literal | carpet_lines_reach(n_max, level)318 if carpet_line_lit(f.numerator, f.denominator, n, level, rs))319 if literal != form:320 bad += 1321 return bad322323def blockmax(a, side):324 f = a.shape[0] // side325 return a.reshape(side, f, side, f).max(axis=(1, 3))326327def grey(a):328 m = int(a.max())329 v = np.zeros(a.shape, dtype=np.float64)330 lit = a > 0331 v[lit] = INK_FLOOR + (255.0 - INK_FLOOR) * (a[lit] - 1) / max(m - 1, 1)332 return (255 - np.round(v)).astype(np.uint8)333334def save_png(acc, path, side):335 Image.fromarray(grey(blockmax(acc, side)), mode="L").save(path, optimize=True)336 return path.stat().st_size, int(acc.max())337338def contact_sheet(accs, path, side):339 tiles = [grey(blockmax(a, side)) for a in accs]340 sheet = np.vstack([np.hstack(tiles[:2]), np.hstack(tiles[2:])])341 Image.fromarray(sheet, mode="L").save(path, optimize=True)342 return path.stat().st_size343344def edge_raster(layers, r):345 acc = np.zeros((r, r), dtype=np.int32)346 for side, v, h in layers:347 for k, j in v:348 c = min(k * r // side, r - 1)349 y0 = j * r // side350 y1 = min(max((j + 1) * r // side, y0 + 1), r)351 acc[y0:y1, c] += 1352 for i, m in h:353 c = min(m * r // side, r - 1)354 x0 = i * r // side355 x1 = min(max((i + 1) * r // side, x0 + 1), r)356 acc[c, x0:x1] += 1357 return acc358359def corner_raster(stack, r, scale):360 acc = np.zeros((r, r), dtype=np.int32)361 for (x, y), b in stack.items():362 cx = min(x.numerator * r // x.denominator, r - 1)363 cy = min(y.numerator * r // y.denominator, r - 1)364 rad = 1 + int(scale * b)365 x0, x1 = max(cx - rad, 0), min(cx + rad + 1, r)366 y0, y1 = max(cy - rad, 0), min(cy + rad + 1, r)367 np.maximum(acc[y0:y1, x0:x1], b, out=acc[y0:y1, x0:x1])368 return acc369370def figures():371 parity_layers = [(n,) + parity_edges(n) for n in range(1, PARITY_RENDER_N + 1, 2)]372 carpet_layers = [(3 ** CARPET_RENDER_L * n,) + carpet_edges(n, CARPET_RENDER_L)373 for n in range(1, CARPET_RENDER_N + 1)]374 accs = [edge_raster(parity_layers, RENDER_R),375 corner_raster(parity_corner_stack(PARITY_RENDER_N), RENDER_R, 1.0),376 edge_raster(carpet_layers, RENDER_R),377 corner_raster(carpet_corner_stack(CARPET_RENDER_N, CARPET_RENDER_L), RENDER_R, 0.25)]378 names = ["edges-parity.png", "corners-parity.png", "edges-carpet.png", "corners-carpet.png"]379 for acc, name in zip(accs, names):380 size, peak = save_png(acc, HERE / name, SAVE_R)381 print("figure", name, "bytes", size, "peak", peak)382 print("figure nodes-sheet.png bytes",383 contact_sheet(accs, HERE / "nodes-sheet.png", SHEET_R))384385def main():386 print("parity edge set literal against form, odd n <= %d, mismatches" % PARITY_CHECK_N,387 check_parity_edges(PARITY_CHECK_N))388 print("parity corner set literal against form, odd n <= %d, mismatches" % PARITY_CHECK_N,389 check_parity_corners(PARITY_CHECK_N))390 literal = parity_corner_stack(PARITY_CHECK_N)391 predicted = parity_corner_predicted(PARITY_CHECK_N)392 print("parity corner stack N = %d, lit points" % PARITY_CHECK_N, len(literal),393 "predicted", len(predicted),394 "missed", len(set(literal) - set(predicted)),395 "invented", len(set(predicted) - set(literal)),396 "value mismatches", sum(1 for p in literal if literal[p] != predicted.get(p)))397 print("parity corner stack digest", stack_digest(literal), stack_digest(predicted))398 top = sorted(literal.items(), key=lambda kv: (-kv[1], kv[0]))[:TOP]399 print("parity corner top", ["%s,%s:%d" % (p[0], p[1], b) for p, b in top])400 print("parity edge lines N = %d, lit lines and form breaches" % PARITY_CHECK_N,401 check_line_brightness_parity(PARITY_CHECK_N))402 print("parity edge segments N = %d, tests and form breaches" % PARITY_CHECK_N,403 parity_segment_check(PARITY_CHECK_N))404 for level in (1, 2):405 js, rows = edge_residues(level)406 print("carpet level %d edge residues J_%d" % (level, level), js,407 "full nonzero", js == list(range(1, 3 ** level)),408 "rows", {j: rows[j] for j in js})409 for n_max, level in CARPET_CHECKS:410 print("carpet edge lines N = %d L = %d, layer mismatches" % (n_max, level),411 check_carpet_lines(n_max, level))412 lit = set()413 for n in range(1, n_max + 1):414 lit |= carpet_edge_lines(n, level)415 reach = carpet_lines_reach(n_max, level)416 print("carpet lit lines N = %d L = %d" % (n_max, level), len(lit),417 "reach form", len(reach), "equal", lit == reach)418 print("carpet line brightness N = %d L = %d, lit lines and form breaches" % (n_max, level),419 check_line_brightness_carpet(n_max, level))420 segs = top_segments(carpet_edge_segments(n_max, level), TOP)421 print("carpet top segments N = %d L = %d" % (n_max, level),422 ["%d at x = %s on [%s, %s]" % (b, x, lo, hi) for (b, w, x, lo, hi) in segs])423 print("carpet edge segments N = %d L = %d, tests and form breaches" % (n_max, level),424 carpet_segment_check(n_max, level))425 stack = carpet_corner_stack(CARPET_CORNER_N, 1)426 bad = sum(1 for p, b in stack.items()427 if carpet_corner_form(p[0].denominator, p[1].denominator, CARPET_CORNER_N) != b)428 print("carpet corner stack N = %d L = 1, lit corners" % CARPET_CORNER_N, len(stack),429 "form mismatches", bad, "digest", stack_digest(stack))430 top = sorted(stack.items(), key=lambda kv: (-kv[1], kv[0]))[:TOP]431 print("carpet corner top", ["%s,%s:%d" % (p[0], p[1], b) for p, b in top])432 tile = carpet_corners(1, 2)433 grid = set((Fraction(k, 9), Fraction(m, 9)) for k in range(10) for m in range(10))434 print("carpet level 2 tile vertices", len(grid), "corners", len(tile),435 "interior to a hole", sorted("%s,%s" % p for p in grid - tile))436 phi = totients(3 * CARPET_CORNER_N)437 print("plain Farey pair count at Q = %d" % (3 * CARPET_CORNER_N),438 (1 + sum(phi[1:])) ** 2)439 print("restricted Farey sizes, L = 1, N = 1..%d" % FAREY_MAX_N,440 [farey_sizes(3 * n) for n in range(1, FAREY_MAX_N + 1)])441 for q in DISCREPANCY_QS:442 plain = landau(farey_list(q))443 restricted = landau(farey_list(q, 3))444 print("landau Q = %d" % q, "plain %.6f" % float(plain),445 "restricted %.6f" % float(restricted))446 figures()447448if __name__ == "__main__":449 main()