"""Part 7: counting the fields within r steps without building them.

The ring construction in hyperbolic.py needs to know one thing about each corner on the boundary:
how many of its five squares are already there. Call the numbers of boundary corners touched by
one, two, three and four squares a, b, c and d. Then the next ring follows four rules:

- a corner touched by one or two squares gets two new fields, stays on the boundary with three or
  four, and sends out two new edges whose far ends are new corners touched by one square;
- a corner touched by three gets two new fields that meet along one new edge, and its far end is a
  new corner touched by two;
- a corner touched by four is a notch: one new field fills it, and the two new edges that would
  have come from its neighbours end at the same corner, so two new corners become one.

So   a' = 2a + 2b - d,   b' = c,   c' = a,   d' = b,
and the next ring has a + b + c fields: one for every corner missing two or more squares. That is
a 4 x 4 matrix, and raising it to the r-th power counts ring r without building anything.
"""
import math
from fractions import Fraction

import hyperbolic

RULE = ((2, 2, 0, -1),
        (0, 0, 1, 0),
        (1, 0, 0, 0),
        (0, 1, 0, 0))
START = (4, 0, 0, 0)                # one square: four corners, each touched by one


def step(v):
    return tuple(sum(RULE[i][j] * v[j] for j in range(4)) for i in range(4))


def mat_mul(a, b):
    return tuple(tuple(sum(a[i][k] * b[k][j] for k in range(4)) for j in range(4)) for i in range(4))


def mat_pow(m, n):
    """m to the n-th power by repeated squaring: about 2 log2(n) multiplications."""
    result = tuple(tuple(int(i == j) for j in range(4)) for i in range(4))
    while n:
        if n & 1:
            result = mat_mul(result, m)
        m = mat_mul(m, m)
        n >>= 1
    return result


def shell(r):
    """Fields exactly r steps from the start, from the matrix alone."""
    if r == 0:
        return 1
    v = mat_pow(RULE, r - 1)
    a, b, c, _ = (sum(v[i][j] * START[j] for j in range(4)) for i in range(4))
    return a + b + c


def ball(r):
    """Fields within r steps: the shells added up."""
    total, v = 1, START
    for _ in range(r):
        total += v[0] + v[1] + v[2]
        v = step(v)
    return total


def boundary_degrees(depth):
    """(a, b, c, d) for the boundary after each ring, read off the actual construction."""
    tiles, layer, _ = hyperbolic.rings(depth)
    out, deg = [], {}
    edges = {}
    for n in range(depth + 1):
        for t in (t for t, l in enumerate(layer) if l == n):
            for k in range(4):
                v = tiles[t][k]
                deg[v] = deg.get(v, 0) + 1
                e = frozenset((v, tiles[t][(k + 1) % 4]))
                edges[e] = edges.get(e, 0) + 1
        on_boundary = {v for e, m in edges.items() if m == 1 for v in e}
        counts = [0, 0, 0, 0]
        for v in on_boundary:
            counts[deg[v] - 1] += 1
        out.append(tuple(counts))
    return out


GROWTH = (1 + math.sqrt(3) + math.sqrt(2 * math.sqrt(3))) / 2


def growth_rate():
    """The largest root of x^4 - 2x^3 - 2x + 1, the matrix's characteristic polynomial. The
    polynomial reads the same backwards, so x + 1/x = y turns it into y^2 - 2y - 2 = 0, and
    y = 1 + sqrt(3) gives x = (1 + sqrt 3 + sqrt(2 sqrt 3)) / 2."""
    return GROWTH


def characteristic():
    """The coefficients of det(xI - RULE), highest power first, by expanding the determinant
    with exact integers (Faddeev-LeVerrier)."""
    n = 4
    ident = tuple(tuple(int(i == j) for j in range(n)) for i in range(n))
    coeffs = [Fraction(1)]
    m = tuple(tuple(0 for _ in range(n)) for _ in range(n))
    for k in range(1, n + 1):
        m = tuple(tuple(sum(RULE[i][t] * m[t][j] for t in range(n)) + int(coeffs[-1]) * ident[i][j]
                        for j in range(n)) for i in range(n))
        am = mat_mul(RULE, m)
        coeffs.append(Fraction(-sum(am[i][i] for i in range(n)), k))
    return [int(c) for c in coeffs]