Fusion rings from groups, and the Fibonacci ring

"""Fusion rings: decategorifying Rep(G), and a fusion ring that is not Rep(G).

The Grothendieck ring of Rep(G) has a basis of irreducible representations,
with V_i (x) V_j = sum_k N_ij^k V_k, and N_ij^k = <chi_i chi_j, chi_k>. We
build S3, D4 and Q8 as groups of matrices, find their one-dimensional
characters by brute force, take traces of a faithful two-dimensional
representation, check that these characters are orthonormal and exhaust the
group (sum of squares of dimensions = |G|), and compute the fusion rings
exactly. D4 and Q8 are not isomorphic (they have 5 and 1 elements of order 2)
but their fusion rings are isomorphic - and Rep(D4), Rep(Q8) are not even
equivalent as tensor categories (Tambara and Yamagami, 1998): the
Grothendieck ring forgets. Finally the Fibonacci fusion ring,
tau (x) tau = 1 + tau, has Frobenius-Perron dimension the golden ratio, so it
is not the Grothendieck ring of any Rep(G); the Ising ring has FPdim sqrt 2.
Exact for the groups; the FP dimensions are floating-point approximations.
"""

from fractions import Fraction
from itertools import permutations, product
import math

# Matrices have Gaussian-integer entries (re, im), so all arithmetic is exact.
ONE = ((1, 0), (0, 0)), ((0, 0), (1, 0))

def g_mul(a, b):
    """Product of Gaussian integers a = (re, im), b = (re, im)."""
    return (a[0] * b[0] - a[1] * b[1], a[0] * b[1] + a[1] * b[0])

def m_mul(x, y):
    """Product of 2 x 2 matrices of Gaussian integers."""
    def entry(i, j):
        p, q = g_mul(x[i][0], y[0][j]), g_mul(x[i][1], y[1][j])
        return (p[0] + q[0], p[1] + q[1])
    return ((entry(0, 0), entry(0, 1)), (entry(1, 0), entry(1, 1)))

def generate(gens):
    """The finite group generated by some invertible 2 x 2 matrices."""
    group, frontier = {ONE}, [ONE]
    while frontier:
        new = []
        for g in frontier:
            for s in gens:
                h = m_mul(g, s)
                if h not in group:
                    group.add(h)
                    new.append(h)
        frontier = new
    return [ONE] + sorted(group - {ONE})  # identity first

def trace(m):
    """Trace, as a Gaussian integer."""
    return (m[0][0][0] + m[1][1][0], m[0][0][1] + m[1][1][1])

def element_order(m):
    """Order of a group element."""
    k, x = 1, m
    while x != ONE:
        x, k = m_mul(x, m), k + 1
    return k

def linear_characters(group):
    """All homomorphisms group -> {1, -1}, by checking every function."""
    index = {g: i for i, g in enumerate(group)}
    table = [[index[m_mul(a, b)] for b in group] for a in group]
    found = []
    for signs in product((1, -1), repeat=len(group)):
        if all(signs[table[i][j]] == signs[i] * signs[j]
               for i in range(len(group)) for j in range(len(group))):
            found.append(list(signs))
    return found

def character_table(gens):
    """Linear characters plus the trace of the defining 2-dim representation.

    For S3, D4 and Q8 the characters here are real integers.
    """
    group = generate(gens)
    chars = linear_characters(group)
    two = [trace(g) for g in group]
    assert all(t[1] == 0 for t in two), "real-valued for these examples"
    chars.append([t[0] for t in two])
    return group, chars

def inner(f, g):
    """<f, g> = average of f * g over the group (real characters)."""
    return Fraction(sum(a * b for a, b in zip(f, g)), len(f))

def fusion(chars):
    """N[i][j][k] = <chi_i chi_j, chi_k>."""
    r = range(len(chars))
    return [[[inner([x * y for x, y in zip(chars[i], chars[j])], chars[k])
              for k in r] for j in r] for i in r]

def check_fusion_ring(n):
    """Non-negative integers, unit = basis element 0, commutative, associative."""
    r = range(len(n))
    integral = all(v.denominator == 1 and v >= 0
                   for a in n for b in a for v in b)
    unit = all(n[0][j][k] == (j == k) for j in r for k in r)
    comm = all(n[i][j][k] == n[j][i][k] for i in r for j in r for k in r)
    assoc = all(sum(n[i][j][m] * n[m][k][l] for m in r)
                == sum(n[j][k][m] * n[i][m][l] for m in r)
                for i, j, k, l in product(r, repeat=4))
    return integral and unit and comm and assoc

def isomorphic_fusion(n1, n2):
    """Is there a relabelling of basis elements (fixing the unit) with n1 = n2?"""
    size = len(n1)
    if size != len(n2):
        return False
    r = range(size)
    for perm in permutations(range(1, size)):
        s = (0,) + perm
        if all(n1[i][j][k] == n2[s[i]][s[j]][s[k]]
               for i in r for j in r for k in r):
            return True
    return False

def fp_dimension(matrix, steps=200):
    """Perron-Frobenius eigenvalue of a non-negative irreducible matrix.

    Power iteration on M + I (whose dominant eigenvalue is strictly largest in
    absolute value, so the iteration converges even when M is periodic, as
    for Ising), then subtract 1.
    """
    size = len(matrix)
    shifted = [[matrix[i][j] + (i == j) for j in range(size)]
               for i in range(size)]
    v = [1.0] * size
    lam = 0.0
    for _ in range(steps):
        w = [sum(a * b for a, b in zip(row, v)) for row in shifted]
        lam = max(w) / max(v)
        v = [x / max(w) for x in w]
    return lam - 1

Z, P, M, I, J = (0, 0), (1, 0), (-1, 0), (0, 1), (0, -1)
GROUPS = {
    # S3 = Weyl group of A2: the two simple reflections, in a basis of roots.
    "S3": [((M, P), (Z, P)), ((P, Z), (P, M))],
    # D4: rotation by 90 degrees and a reflection.
    "D4": [((Z, M), (P, Z)), ((P, Z), (Z, M))],
    # Q8: the quaternions i and j as 2 x 2 complex matrices.
    "Q8": [((I, Z), (Z, J)), ((Z, P), (M, Z))],
}

def main():
    rings = {}
    for name, gens in GROUPS.items():
        group, chars = character_table(gens)
        n = fusion(chars)
        rings[name] = n
        dims = [c[0] for c in chars]
        orth = all(inner(a, b) == (i == j) for i, a in enumerate(chars)
                   for j, b in enumerate(chars))
        invol = sum(1 for g in group if element_order(g) == 2)
        squares = sum(d * d for d in dims)
        ok = orth and check_fusion_ring(n) and squares == len(group)
        print(f"Rep({name}): |G| = {len(group)}, dims {dims}, elements of "
              f"order 2: {invol}; checks {'pass' if ok else 'FAIL'}")
    last = len(rings["S3"]) - 1
    print("  in Rep(S3): V (x) V has multiplicities",
          [int(v) for v in rings["S3"][last][last]], "(triv, sign, V)")
    print("  in Rep(D4): E (x) E has multiplicities",
          [int(v) for v in rings["D4"][4][4]], "(four lines, E)")
    print("  fusion rings of D4 and Q8 isomorphic:",
          isomorphic_fusion(rings["D4"], rings["Q8"]))
    print()
    fib = [[0, 1], [1, 1]]  # multiplication by tau in the basis (1, tau)
    lam = fp_dimension(fib)
    print("Fibonacci ring: tau (x) tau = 1 + tau")
    print(f"  FPdim(tau) by power iteration = {lam:.12f}")
    print(f"  golden ratio (1 + sqrt 5)/2    = {(1 + math.sqrt(5)) / 2:.12f}")
    ising = [[0, 1, 0], [1, 0, 1], [0, 1, 0]]  # sigma on (1, sigma, psi)
    print("Ising ring: sigma (x) sigma = 1 + psi, sigma (x) psi = sigma")
    print(f"  FPdim(sigma) = {fp_dimension(ising):.12f}"
          f" (sqrt 2 = {math.sqrt(2):.12f})")
    print()
    print("Epistemic status: the groups, characters, fusion coefficients and")
    print("ring axioms are computed exactly and exhaustively for these three")
    print("groups. The FP dimension is a floating-point approximation; that")
    print("d^2 = 1 + d has no integer root is proved in the lesson, and that")
    print("Rep(D4), Rep(Q8) differ as tensor categories is a theorem not")
    print("checked here.")

if __name__ == "__main__":
    main()