"""Batch 3: verify the best survivors much further, add remaining areas."""
from math import isqrt, gcd
from itertools import combinations


def say(s=""):
    print(s)


def det_bareiss(M):
    M = [r[:] for r in M]
    n = len(M)
    if n == 0:
        return 1
    sign, prev = 1, 1
    for k in range(n - 1):
        if M[k][k] == 0:
            piv = next((r for r in range(k + 1, n) if M[r][k] != 0), None)
            if piv is None:
                return 0
            M[k], M[piv] = M[piv], M[k]
            sign = -sign
        for i in range(k + 1, n):
            for j in range(k + 1, n):
                M[i][j] = (M[i][j] * M[k][k] - M[i][k] * M[k][j]) // prev
        prev = M[k][k]
    return sign * M[n - 1][n - 1]


# ---------- A. Hankel determinants, deeper and wider ----------
L = 400


def rudin_shapiro(n):
    return (-1) ** sum(1 for k in range(24) if (n >> k) & 3 == 3)


def baum_sweet(n):
    if n == 0:
        return 1
    b = bin(n)[2:]
    run = 0
    for ch in b:
        if ch == "0":
            run += 1
        else:
            if run % 2 == 1:
                return 0
            run = 0
    return 0 if run % 2 == 1 else 1


def stern(n):
    a, b = 1, 0
    while n:
        if n & 1:
            b += a
        else:
            a += b
        n >>= 1
    return b


def period_doubling(n):
    m = n + 1
    return 1 if (m // (m & -m)) % 4 == 1 else 0


def cantor(n):
    while n:
        if n % 3 == 1:
            return 0
        n //= 3
    return 1


SEQ = {
    "rudin_shapiro": [rudin_shapiro(n) for n in range(L)],
    "baum_sweet": [baum_sweet(n) for n in range(L)],
    "stern": [stern(n) for n in range(L)],
    "period_doubling": [period_doubling(n) for n in range(L)],
    "cantor_ternary": [cantor(n) for n in range(L)],
    "thue_morse_pm1": [1 - 2 * (bin(n).count("1") % 2) for n in range(L)],
}

say("=== A. Hankel determinants det(a_{i+j})_{0<=i,j<n}, n up to 80 ===")
say(f"{'sequence':18} {'zeros (n<=80)':>40}")
say("-" * 62)
for name, a in SEQ.items():
    zeros = []
    for n in range(1, 81):
        if det_bareiss([[a[i + j] for j in range(n)] for i in range(n)]) == 0:
            zeros.append(n)
    say(f"{name:18} {str(zeros[:14]):>40}")

# ---------- B. permutations with i + sigma(i) always a square / cube ----------
say()
say("=== B. Permutations sigma of {1..n} with i+sigma(i) always a perfect square/cube ===")


def perfect_matching_exists(n, pred):
    adj = {i: [j for j in range(1, n + 1) if pred(i + j)] for i in range(1, n + 1)}
    matchR = [-1] * (n + 1)
    def aug(u, seen):
        for v in adj[u]:
            if not seen[v]:
                seen[v] = True
                if matchR[v] == -1 or aug(matchR[v], seen):
                    matchR[v] = u
                    return True
        return False
    return all(aug(u, [False] * (n + 1)) for u in range(1, n + 1))


sq = lambda m: isqrt(m) ** 2 == m
cu = lambda m: round(m ** (1 / 3)) ** 3 == m
tri = lambda m: isqrt(8 * m + 1) ** 2 == 8 * m + 1
for label, pred in (("square", sq), ("cube", cu), ("triangular", tri)):
    ok = [n for n in range(1, 81) if perfect_matching_exists(n, pred)]
    bad = [n for n in range(1, 81) if n not in ok]
    say(f"  i+sigma(i) always a {label:11}: works for {len(ok)}/80; fails at n = {bad[:18]}")

# ---------- C. domino tilings of 4 x n mod small primes ----------
say()
say("=== C. Domino tilings of a 4 x n rectangle, residues mod small primes ===")
T = [1, 1, 5, 11, 36, 95, 281, 781, 2245, 6336, 18061, 51205, 145601, 413351,
     1174500, 3335651, 9475901, 26915305, 76455961, 217172736, 616891945]
for p in (2, 3, 5, 7):
    r = [t % p for t in T]
    per = None
    for q in range(1, 11):
        if all(r[k] == r[k + q] for k in range(len(r) - q)):
            per = q
            break
    say(f"  mod {p}: {r}   period<=10: {per}")

# ---------- D. grid points: no 3 collinear and no 4 concyclic ----------
say()
say("=== D. Max points in n x n grid, no 3 collinear and no 4 concyclic ===")


def collinear(p, q, r):
    return (q[0] - p[0]) * (r[1] - p[1]) == (q[1] - p[1]) * (r[0] - p[0])


def concyc(p, q, r, s):
    def d3(a):
        return (a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
                - a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0])
                + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]))
    M = [[x * x + y * y, x, y, 1] for (x, y) in (p, q, r, s)]
    tot = 0
    for c in range(4):
        minor = [[M[r2][c2] for c2 in range(4) if c2 != c] for r2 in range(1, 4)]
        tot += ((-1) ** c) * M[0][c] * d3(minor)
    return tot == 0


def best_grid(n):
    pts = [(x, y) for x in range(n) for y in range(n)]
    best = 0
    chosen = []
    def ok(c):
        for a, b in combinations(chosen, 2):
            if collinear(a, b, c):
                return False
        for trip in combinations(chosen, 3):
            if concyc(*trip, c):
                return False
        return True
    def bt(idx):
        nonlocal best
        best = max(best, len(chosen))
        if idx == len(pts) or len(chosen) + (len(pts) - idx) <= best:
            return
        c = pts[idx]
        if ok(c):
            chosen.append(c)
            bt(idx + 1)
            chosen.pop()
        bt(idx + 1)
    bt(0)
    return best


for n in range(2, 8):
    say(f"  n={n}: {best_grid(n)}")