#!/usr/bin/env python3
"""Bounded exact checks for the relative and cyclic Postnikov examples.

Python standard library only; no external packages, randomness, or floats.
Run: python3 verify_relative_examples.py > relative_examples_results.json

Bounds: normalized-bar matrices over F_p for (p,r)=(3,2),(3,3),(3,4),
(5,2), with every (a,b,c_Y,c_X) in F_p^4; twisted C_2 cochains over
Z/4 through degree 10 and lift parity for 2 <= r <= 8; odd-prime
binomial identities for p=3,5,7, 2 <= j <= 100, every nonzero a;
central-extension multiplication tables for p=3,5,7 and every m.

These are finite consistency checks, NOT proofs in arbitrary degree or
verification of the full homotopy types. Nontriviality and solvability
are computed from actual normalized-bar coboundary images by modular
Gaussian elimination, independently of the scalar criterion tested.
"""

from itertools import product
from math import comb
import hashlib
import json
from pathlib import Path


def tuples(p, degree):
    return list(product(range(1, p), repeat=degree))


def differential_terms(p, args, action=lambda g: 1):
    """Signed inhomogeneous bar differential, before normalization."""
    n = len(args) - 1
    yield action(args[0]), args[1:]
    for i in range(1, n + 1):
        merged = args[:i-1] + ((args[i-1] + args[i]) % p,) + args[i+1:]
        yield (-1) ** i, merged
    yield (-1) ** (n + 1), args[:-1]


def bar_matrix(p, degree):
    source = tuples(p, degree)
    target = tuples(p, degree + 1)
    index = {args: i for i, args in enumerate(source)}
    rows = []
    for args in target:
        row = [0] * len(source)
        for sign, face in differential_terms(p, args):
            if 0 not in face:
                col = index[face]
                row[col] = (row[col] + sign) % p
        rows.append(row)
    return rows


def matvec(matrix, vector, p):
    return [sum(x*y for x, y in zip(row, vector)) % p for row in matrix]


def column_basis(matrix, p):
    """Echelon basis for the image, stored as (pivot, normalized vector)."""
    basis = []
    for col in zip(*matrix):
        v = list(col)
        for pivot, b in basis:
            factor = v[pivot]
            if factor:
                v = [(x-factor*y) % p for x, y in zip(v, b)]
        pivot = next((i for i, x in enumerate(v) if x), None)
        if pivot is not None:
            inv = pow(v[pivot], -1, p)
            basis.append((pivot, [(x*inv) % p for x in v]))
            basis.sort(key=lambda pair: pair[0])
    return basis


def in_image(vector, basis, p):
    v = list(vector)
    for pivot, b in basis:
        factor = v[pivot]
        if factor:
            v = [(x-factor*y) % p for x, y in zip(v, b)]
    return not any(v)


def xi(p, args):
    """Actual cup cochain v cup u^j (odd) or u^j (even)."""
    if 0 in args:
        return 0
    start = len(args) % 2
    value = args[0] if start else 1
    for i in range(start, len(args), 2):
        value *= (args[i] + args[i+1]) // p
    return value % p


def cyclic_case(p, r):
    d = r + 1
    matrices = {q: bar_matrix(p, q) for q in range(d + 1)}
    bases = {q: column_basis(matrices[q], p) for q in matrices}
    cohomology = []
    for q in range(d + 1):
        previous_rank = len(bases[q-1]) if q else 0
        rank = len(bases[q])
        dimension = (p-1)**q - previous_rank - rank
        assert dimension == 1
        cocycle = [xi(p, args) for args in tuples(p, q)]
        assert not any(matvec(matrices[q], cocycle, p))
        assert not in_image(cocycle, bases[q-1] if q else [], p)
        # Check d^2=0 on every normalized basis cochain in the range.
        if q:
            for col in zip(*matrices[q-1]):
                assert not any(matvec(matrices[q], col, p))
        cohomology.append({
            'degree': q, 'cochain_dimension': (p-1)**q,
            'incoming_rank': previous_rank, 'outgoing_rank': rank,
            'cohomology_dimension': dimension,
            'xi_is_cocycle_and_not_boundary': True,
        })

    args_d = tuples(p, d)
    invariant = [xi(p, args) for args in args_d]
    pullbacks = {a: [xi(p, tuple(a*g % p for g in args))
                     for args in args_d] for a in range(p)}
    e = (r + 2) // 2
    for a in range(p):
        residual = [(x - pow(a, e, p)*y) % p
                    for x, y in zip(pullbacks[a], invariant)]
        assert in_image(residual, bases[d-1], p)

    soluble = 0
    for a, b, cy, cx in product(range(p), repeat=4):
        obstruction = [(cx*x-b*cy*y) % p
                       for x, y in zip(pullbacks[a], invariant)]
        bar_answer = in_image(obstruction, bases[d-1], p)
        predicted = (pow(a, e, p)*cx-b*cy) % p == 0
        assert bar_answer == predicted, (p, r, a, b, cy, cx)
        soluble += bar_answer
    # Integral lift of v establishes the sign of its mod-p Bockstein.
    for g, h in product(range(p), repeat=2):
        dv = h - ((g+h) % p) + g
        assert dv == p*((g+h)//p)
    return {
        'p': p, 'r': r, 'obstruction_degree': d,
        'all_scalar_quadruples_checked': p**4,
        'solvable_quadruples': soluble,
        'scalar_equation_matches_bar_image': True,
        'positive_carry_bockstein_sign_checked': True,
        'cohomology_by_actual_bar_ranks': cohomology,
    }


def twisted_parity():
    # There is one normalized C2 tuple in each degree. Evaluate the
    # actual bar formula with the sign action, not an imposed parity rule.
    rows = []
    maps = {}
    for degree in range(11):
        args = (1,)*(degree+1)
        values = []
        for value in range(4):
            delta = sum(sign*value for sign, face in
                        differential_terms(2, args, lambda g: (-1)**g)
                        if 0 not in face) % 4
            values.append(delta)
        maps[degree] = values
        assert values == [(2 if degree % 2 == 0 else 0)*a % 4
                          for a in range(4)]
    for degree in range(10):
        kernel = {x for x in range(4) if maps[degree][x] == 0}
        image = set(maps[degree-1]) if degree else {0}
        assert image <= kernel
        assert len(kernel)//len(image) == 2
        representatives = []
        remaining = set(kernel)
        while remaining:
            rep = min(remaining)
            representatives.append(rep)
            remaining -= {(rep+x) % 4 for x in image}
        nonzero = next(rep for rep in representatives if rep not in image)
        reduction = nonzero % 2
        assert reduction == degree % 2
        rows.append({'degree': degree, 'kernel': sorted(kernel),
                     'boundary_image': sorted(image),
                     'nonzero_class_representative': nonzero,
                     'nonzero_class_mod_2': reduction})
    lifts = []
    for r in range(2, 9):
        obstruction = rows[r+1]['nonzero_class_mod_2']
        assert (obstruction == 0) == (r % 2 == 1)
        lifts.append({'r': r, 'obstruction_mod_2': obstruction,
                      'lift_exists': obstruction == 0})
    return {'group': 'C2', 'band': 'Z/4 with sign action',
            'differential_degree_max': 10, 'cohomology': rows,
            'lift_parity_r_min': 2, 'lift_parity_r_max': 8,
            'lifts': lifts}


def padic_coefficients():
    cases = 0
    boundary_examples = []
    for p in (3, 5, 7):
        for j in range(2, 101):
            q, s = j, 0
            while q % p == 0:
                q //= p
                s += 1
            first = next(l for l in range(1, j+1) if comb(j, l) % p)
            assert first == p**s
            for a in range(1, p):
                coefficients = {l: comb(j, l)*pow(a, j-l, p) % p
                                for l in range(1, j+1)}
                for l, coefficient in coefficients.items():
                    if l % (p**s):
                        assert coefficient == 0
                    else:
                        expected = comb(q, l//(p**s))*pow(a, j-l, p) % p
                        assert coefficient == expected
                assert coefficients[first] == q*pow(a, j-first, p) % p
                cases += 1
            if j in (p, p*p, 2*p):
                boundary_examples.append({'p': p, 'j': j, 'valuation': s,
                                          'first_nonzero_exponent': first,
                                          'retained_homotopy_degree': 2*first-1})
    return {'primes': [3, 5, 7], 'j_min': 2, 'j_max': 100,
            'all_nonzero_a_checked': True, 'number_of_p_j_a_cases': cases,
            'all_coefficient_identities_hold': True,
            'boundary_examples': boundary_examples}


def central_extensions():
    rows = []
    for p in (3, 5, 7):
        elements = list(product(range(p), repeat=2))
        for m in range(p):
            def mul(x, y):
                return ((x[0]+y[0]) % p,
                        (x[1]+y[1]+m*((x[0]+y[0])//p)) % p)
            for x, y in product(elements, repeat=2):
                assert mul(x, y) == mul(y, x)
            for x, y, z in product(elements, repeat=3):
                assert mul(mul(x, y), z) == mul(x, mul(y, z))
            value, order = (0, 0), 0
            while True:
                value = mul(value, (1, 0))
                order += 1
                if value == (0, 0):
                    break
            assert order == (p if m == 0 else p*p)
            rows.append({'p': p, 'm': m, 'order_of_1_0': order,
                         'group': 'Cp x Cp' if m == 0 else 'C(p^2)'})
    return {'primes': [3, 5, 7], 'every_m_checked': True,
            'full_associativity_tables_checked': True, 'cases': rows}


def main():
    report = {
        'method': 'Exact normalized bar matrices and modular Gaussian elimination; integer enumeration.',
        'scope_warning': 'Finite consistency checks only. General theorems and space-level splittings require the manuscript proofs.',
        'cyclic_lifts': [cyclic_case(p, r) for p, r in ((3, 2), (3, 3), (3, 4), (5, 2))],
        'twisted_C2_parity': twisted_parity(),
        'padic_gauge_coefficients': padic_coefficients(),
        'fundamental_group_extensions': central_extensions(),
        'script_sha256': hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
    }
    print(json.dumps(report, indent=2, sort_keys=True))


if __name__ == '__main__':
    main()