Quotient groups and quotient rings

"""When is "multiply the cosets" a well-defined operation?

For every subgroup H of the groups S3, D4 and A4 (as permutation groups), the
program tests the rule (aH)(bH) = abH against every choice of representatives,
and separately tests whether H is normal (gHg^-1 = H for all g). The two tests
always agree, as the lesson proves. For each normal subgroup it reports the
order of the quotient group G/H and the orders of its elements; for each
non-normal one it prints a witness showing the rule depends on the choice of
representative. Exact and exhaustive for these three groups.
"""

from collections import Counter

def compose(p, q):
    """The permutation 'first q, then p'."""
    return tuple(p[i] for i in q)

def inverse(p):
    """Inverse permutation."""
    inv = [0] * len(p)
    for i, image in enumerate(p):
        inv[image] = i
    return tuple(inv)

def closure(gens):
    """The group generated by a list of permutations of the same degree."""
    identity = tuple(range(len(gens[0])))
    elements, frontier = {identity}, [identity]
    while frontier:
        new = []
        for x in frontier:
            for g in gens:
                y = compose(g, x)
                if y not in elements:
                    elements.add(y)
                    new.append(y)
        frontier = new
    return frozenset(elements)

def subgroups(group):
    """All subgroups, as iterated joins of cyclic subgroups (exhaustive)."""
    found = {closure([g]) for g in group}
    frontier = list(found)
    while frontier:
        new = []
        for h in frontier:
            for k in list(found):
                join = closure(sorted(h | k))
                if join not in found:
                    found.add(join)
                    new.append(join)
        frontier = new
    return sorted(found, key=lambda s: (len(s), sorted(s)))

def coset(a, sub):
    """The left coset aH as a frozenset."""
    return frozenset(compose(a, h) for h in sub)

def is_normal(sub, group):
    """gHg^-1 = H for every g in the group."""
    return all(compose(compose(g, h), inverse(g)) in sub
               for g in group for h in sub)

def product_witness(sub, group):
    """Return (a, a2, b) with aH = a2H but abH != a2bH, or None if the
    coset product is well defined. Left factors are varied; right factors
    are fixed, which already detects every failure (see the lesson)."""
    for a in sorted(group):
        for h in sorted(sub):
            a2 = compose(a, h)
            for b in sorted(group):
                if coset(compose(a, b), sub) != coset(compose(a2, b), sub):
                    return a, a2, b
    return None

def quotient_orders(sub, group):
    """Orders of the elements of G/N, computed on cosets."""
    cosets = {coset(g, sub) for g in group}
    identity = frozenset(sub)
    orders = []
    for c in cosets:
        rep = next(iter(c))
        power, k = rep, 1
        while coset(power, sub) != identity:
            power, k = compose(rep, power), k + 1
        orders.append(k)
    return sorted(Counter(orders).items())

def groups():
    """The three groups, as (name, generators)."""
    return [
        ("S3", [(1, 0, 2), (1, 2, 0)]),
        ("D4", [(1, 2, 3, 0), (0, 3, 2, 1)]),
        ("A4", [(1, 2, 0, 3), (1, 0, 3, 2)]),
    ]

def main():
    agree = True
    witness_shown = False
    for name, gens in groups():
        group = closure(gens)
        print(f"{name} (order {len(group)}): |H|  normal  product defined  G/H")
        for sub in subgroups(group):
            normal = is_normal(sub, group)
            witness = product_witness(sub, group)
            agree &= normal == (witness is None)
            info = f"order {len(group) // len(sub)}, orders "
            info += str(quotient_orders(sub, group)) if normal else "-"
            print(f"   {len(sub):>3}  {str(normal):<6}  {str(witness is None):<15}"
                  f"  {info if normal else 'not a group'}")
            if witness and not witness_shown:
                witness_shown = True
                a, a2, b = witness
                print(f"     e.g. a={a}, a'={a2} lie in one coset, but ab and")
                print(f"     a'b lie in different cosets (b={b})")
    print(f"'product well defined' agrees with 'normal' every time: {agree}")
    print()
    print("Epistemic status: exact and exhaustive for every subgroup of S3, D4")
    print("and A4. It illustrates the theorem (well defined <=> normal), which")
    print("the lesson proves for all groups; the witnesses are genuine proofs")
    print("that those particular subgroups are not normal.")

if __name__ == "__main__":
    main()

"""Building finite fields as quotients F_p[x]/(f), and watching some fail.

Elements of F_p[x]/(f), for f of degree n, are the p^n polynomials of degree
below n, multiplied by reducing modulo f. The program checks, exhaustively,
whether every non-zero element has an inverse (so the quotient is a field) and
whether f has a root or factors, and confirms in each case: field exactly when f
is irreducible. For the fields it finds an element whose powers give every
non-zero element (the multiplicative group is cyclic), and it checks the
Chinese remainder isomorphism Z/12 = Z/3 x Z/4 against Z/4 versus Z/2 x Z/2.
"""

from itertools import product

def poly_mulmod(a, b, f, p):
    """Multiply coefficient lists a, b (lowest degree first) modulo f and p.

    f is monic of degree n; the result has length n."""
    n = len(f) - 1
    prod = [0] * (2 * n)
    for i, x in enumerate(a):
        for j, y in enumerate(b):
            prod[i + j] = (prod[i + j] + x * y) % p
    for k in range(2 * n - 1, n - 1, -1):      # reduce using x^n = -(lower terms)
        c = prod[k]
        if c:
            for t in range(n + 1):
                prod[k - n + t] = (prod[k - n + t] - c * f[t]) % p
    return tuple(prod[:n])

def elements(p, n):
    """All p^n residues: coefficient tuples of length n."""
    return [tuple(c) for c in product(range(p), repeat=n)]

def units(f, p):
    """Residues with a multiplicative inverse modulo f."""
    n = len(f) - 1
    one = tuple([1] + [0] * (n - 1))
    elems = elements(p, n)
    return [a for a in elems if any(poly_mulmod(a, b, f, p) == one for b in elems)]

def has_root(f, p):
    """True when f has a root in F_p."""
    return any(sum(c * pow(x, k, p) for k, c in enumerate(f)) % p == 0
               for x in range(p))

def is_irreducible(f, p):
    """Exhaustive: f (monic, degree n) has no monic factor of degree 1..n//2."""
    n = len(f) - 1
    for d in range(1, n // 2 + 1):
        for tail in product(range(p), repeat=d):
            g = list(tail) + [1]
            # long division of f by g over F_p
            r = list(f)
            for k in range(n - d, -1, -1):
                c = r[k + d]
                for t in range(d + 1):
                    r[k + t] = (r[k + t] - c * g[t]) % p
            if not any(r[:d]):
                return False
    return True

def generator(f, p):
    """A residue whose powers run through all non-zero residues, or None."""
    n = len(f) - 1
    one = tuple([1] + [0] * (n - 1))
    size = p ** n - 1
    for g in elements(p, n):
        if not any(g):
            continue
        power, k = g, 1
        while power != one and k <= size:
            power, k = poly_mulmod(power, g, f, p), k + 1
        if k == size:
            return g
    return None

def show_poly(f):
    """Coefficients listed lowest degree first, as text: [1, 1, 0, 1] -> x^3+x+1."""
    terms = []
    for k in range(len(f) - 1, -1, -1):
        c = f[k]
        if not c:
            continue
        mono = "" if k == 0 else ("x" if k == 1 else f"x^{k}")
        coef = str(c) if (c != 1 or k == 0) else ""
        terms.append(coef + mono)
    return "+".join(terms) if terms else "0"

def crt_is_isomorphism(m, n):
    """x -> (x mod m, x mod n) is a bijective ring map Z/mn -> Z/m x Z/n."""
    images = {(x % m, x % n) for x in range(m * n)}
    mult = all(((x * y) % m, (x * y) % n) == ((x % m) * (y % m) % m,
                                              (x % n) * (y % n) % n)
               for x in range(m * n) for y in range(m * n))
    return len(images) == m * n and mult

CASES = [(2, [1, 1, 1]), (2, [1, 0, 1]), (3, [1, 0, 1]), (3, [1, 1, 1]),
         (2, [1, 1, 0, 1]), (2, [1, 1, 0, 0, 1]), (5, [2, 0, 1])]

def main():
    print(" p  f            irreducible  root  units / non-zero  field  generator")
    for p, f in CASES:
        n = len(f) - 1
        u = units(f, p)
        field = len(u) == p ** n - 1
        g = generator(f, p) if field else None
        print(f"{p:>2}  {show_poly(f):<12} {str(is_irreducible(f, p)):<12} "
              f"{str(has_root(f, p)):<5} {len(u):>5} / {p ** n - 1:<9} "
              f"{str(field):<6} {show_poly(g) if g else '-'}")
    print("Field <=> f irreducible in every case:",
          all((len(units(f, p)) == p ** (len(f) - 1) - 1) == is_irreducible(f, p)
              for p, f in CASES))
    print(f"CRT: Z/12 -> Z/3 x Z/4 is an isomorphism: {crt_is_isomorphism(3, 4)}")
    print(f"     Z/4 -> Z/2 x Z/2 is an isomorphism: {crt_is_isomorphism(2, 2)}")
    print()
    print("Epistemic status: exact and exhaustive for the seven quotients and")
    print("the two maps shown. That k[x]/(f) is a field exactly when f is")
    print("irreducible, and that the unit group of a finite field is cyclic,")
    print("are theorems; these computations illustrate them, they do not prove")
    print("them.")

if __name__ == "__main__":
    main()