"""Batch 2: more areas, plus a deeper look at the survivors from batch 1."""
from math import isqrt
from itertools import combinations

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


def is_prime(m):
    if m < 2:
        return False
    if m % 2 == 0:
        return m == 2
    for p in range(3, isqrt(m) + 1, 2):
        if m % p == 0:
            return False
    return True


# ============ 5. Hankel determinants of automatic sequences ============
# H_n = det( a_{i+j} )_{0<=i,j<n}. For many automatic sequences these are
# conjecturally never zero, which is the same flavour of question as the
# Fibonacci determinant but a different object.
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]


L = 200
thue = [bin(n).count("1") % 2 for n in range(L)]
rudin = [(-1) ** sum(1 for k in range(20) if (n >> k) & 3 == 3) for n in range(L)]


def paperfold(n):
    n += 1
    while n % 2 == 0:
        n //= 2
    return 1 if n % 4 == 1 else 0


paper = [paperfold(n) for n in range(L)]
period_doubling = [1 if (lambda m: (m // (m & -m)) % 4 == 1)(n + 1) else 0 for n in range(L)]
seqs = {
    "thue_morse_01": thue,
    "thue_morse_pm1": [1 - 2 * t for t in thue],
    "rudin_shapiro": rudin,
    "paperfolding": paper,
    "period_doubling": period_doubling,
}
say("=== 5. Hankel determinants det(a_{i+j}) of automatic sequences ===")
say(f"{'sequence':20} {'n with H_n=0 (n<=45)':>26}  max|H_n|")
say("-" * 70)
for name, a in seqs.items():
    zeros, mx = [], 0
    for n in range(1, 46):
        H = det_bareiss([[a[i + j] for j in range(n)] for i in range(n)])
        if H == 0:
            zeros.append(n)
        mx = max(mx, abs(H))
    say(f"{name:20} {str(zeros[:9]):>26}  {mx}")

# ============ 6. Multiplicative Sidon: distinct pairwise products ============
say()
say("=== 6. Largest subset of {1..n} with all pairwise products distinct ===")


def max_mult_sidon(n):
    best = 0
    chosen = []
    def ok(c):
        prods = set()
        for a, b in combinations(chosen + [c], 2):
            p = a * b
            if p in prods:
                return False
            prods.add(p)
        return True
    def bt(v):
        nonlocal best
        best = max(best, len(chosen))
        if v > n or len(chosen) + (n - v + 1) <= best:
            return
        if ok(v):
            chosen.append(v)
            bt(v + 1)
            chosen.pop()
        bt(v + 1)
    bt(1)
    return best


say("  " + "  ".join(f"{n}:{max_mult_sidon(n)}" for n in range(2, 26)))

# ============ 7. Binary words avoiding abelian squares ============
say()
say("=== 7. Binary words with no two adjacent blocks that are anagrams (abelian squares) ===")


def count_abelian_square_free(n):
    def has_ab_sq(w):
        L = len(w)
        for l in range(1, L // 2 + 1):
            a = w[L - 2 * l:L - l]
            b = w[L - l:]
            if sorted(a) == sorted(b):
                return True
        return False
    counts = [0] * (n + 1)
    def rec(w):
        counts[len(w)] += 1
        if len(w) == n:
            return
        for c in (0, 1):
            nw = w + (c,)
            if not has_ab_sq(nw):
                rec(nw)
    rec(())
    return counts


c = count_abelian_square_free(18)
say("  length : count of abelian-square-free binary words")
say("  " + "  ".join(f"{i}:{c[i]}" for i in range(1, 19)))

# ============ 8. Cop number / burning of circulant graphs ============
say()
say("=== 8. Burning number of the circulant graph C_n(1,k) ===")


def burning_number(n, conns):
    adj = [[] for _ in range(n)]
    for v in range(n):
        for d in conns:
            adj[v].append((v + d) % n)
            adj[v].append((v - d) % n)
    from collections import deque
    def dists(src):
        D = [-1] * n
        D[src] = 0
        q = deque([src])
        while q:
            u = q.popleft()
            for w in adj[u]:
                if D[w] < 0:
                    D[w] = D[u] + 1
                    q.append(w)
        return D
    DIST = [dists(v) for v in range(n)]
    for b in range(1, 8):
        found = False
        for centers in combinations(range(n), b):
            cov = set()
            for idx, cvert in enumerate(centers):
                r = b - 1 - idx
                cov |= {u for u in range(n) if DIST[cvert][u] <= r}
            if len(cov) == n:
                found = True
                break
        if found:
            return b
    return None


for k in (2, 3):
    vals = [(n, burning_number(n, (1, k))) for n in range(6, 19)]
    say(f"  C_n(1,{k}): " + "  ".join(f"{n}:{v}" for n, v in vals))

# ============ 9. Deeper look at surviving subtraction games ============
say()
say("=== 9. Subtraction games: does the P-position density go to zero? ===")


def grundy_P(S, N):
    g = [0] * (N + 1)
    P = []
    for n in range(1, N + 1):
        seen = set()
        for m in S:
            if m > n:
                break
            seen.add(g[n - m])
        v = 0
        while v in seen:
            v += 1
        g[n] = v
        if v == 0:
            P.append(n)
    return g, P


N2 = 200000
tri = [m for m in range(1, N2 + 1) if isqrt(8 * m + 1) ** 2 == 8 * m + 1]
cub = [m for m in range(1, N2 + 1) if round(m ** (1 / 3)) ** 3 == m]
sqs = [m for m in range(1, N2 + 1) if isqrt(m) ** 2 == m]
for name, S in (("triangular", tri), ("cubes", cub), ("squares", sqs)):
    g, P = grundy_P(S, N2)
    say(f"  subtract-a-{name}: max grundy {max(g)}, |P| = {len(P)}")
    for cut in (10000, 50000, 100000, 200000):
        d = sum(1 for p in P if p <= cut) / cut
        say(f"      density up to {cut:>7}: {d:.5f}")