"""Exact small calculations used only to assess proposed research projects.

These examples and invariants are classical.  This script makes no novelty claim.
"""

import itertools
import json
import math
from pathlib import Path

import sympy as sp
from sympy.matrices.normalforms import smith_normal_form


def divides(a, b):
    return all(x <= y for x, y in zip(a, b))


def monomial_hom_dimension(generators):
    n = len(generators[0])
    bounds = [
        min(g[i] for g in generators if g[i] and sum(g) == g[i])
        for i in range(n)
    ]
    basis = [
        b
        for b in itertools.product(*(range(c) for c in bounds))
        if not any(divides(g, b) for g in generators)
    ]
    index = {b: i for i, b in enumerate(basis)}
    width = len(basis)
    rows = []
    for i, j in itertools.combinations(range(len(generators)), 2):
        lcm = tuple(max(a, b) for a, b in zip(generators[i], generators[j]))
        output_rows = {}
        for source, sign in [(i, 1), (j, -1)]:
            multiplier = tuple(a - b for a, b in zip(lcm, generators[source]))
            for k, b in enumerate(basis):
                product = tuple(a + c for a, c in zip(multiplier, b))
                if product in index:
                    row = output_rows.setdefault(index[product], {})
                    column = source * width + k
                    row[column] = row.get(column, 0) + sign
        rows.extend(row for row in output_rows.values() if any(row.values()))
    variables = len(generators) * width
    matrix = sp.zeros(len(rows), variables)
    for r, row in enumerate(rows):
        for column, coefficient in row.items():
            matrix[r, column] = coefficient
    hilbert_function = [sum(sum(b) == i for b in basis) for i in range(max(map(sum, basis)) + 1)]
    return {
        "standard_monomial_exponents": basis,
        "length": width,
        "hilbert_function": hilbert_function,
        "hom_dimension_over_Q": variables - matrix.rank(),
    }


def rank_mod_two(matrix):
    rows = [list(map(lambda x: int(x) % 2, row)) for row in matrix.tolist()]
    rank = 0
    for c in range(matrix.cols):
        pivot = next((r for r in range(rank, len(rows)) if rows[r][c]), None)
        if pivot is None:
            continue
        rows[pivot], rows[rank] = rows[rank], rows[pivot]
        for r in range(rank + 1, len(rows)):
            if rows[r][c]:
                rows[r] = [x ^ y for x, y in zip(rows[r], rows[rank])]
        rank += 1
    return rank


def projective_plane_example():
    facets = [
        (1, 2, 3), (1, 2, 4), (1, 3, 5), (1, 4, 6), (1, 5, 6),
        (2, 3, 6), (2, 4, 5), (2, 5, 6), (3, 4, 5), (3, 4, 6),
    ]
    faces = {r: sorted({face for f in facets for face in itertools.combinations(f, r + 1)}) for r in range(3)}
    boundary = {}
    for d in [1, 2]:
        lower_index = {f: i for i, f in enumerate(faces[d - 1])}
        matrix = sp.zeros(len(faces[d - 1]), len(faces[d]))
        for column, f in enumerate(faces[d]):
            for i in range(d + 1):
                matrix[lower_index[f[:i] + f[i + 1:]], column] = (-1) ** i
        boundary[d] = matrix
    zero_product = boundary[1] * boundary[2] == sp.zeros(6, 10)
    output = {"facets": facets, "boundary_squared_zero": zero_product}
    for label, rank_fn in [("Q", lambda m: m.rank()), ("F2", rank_mod_two)]:
        r1, r2 = (rank_fn(boundary[d]) for d in [1, 2])
        output[label] = {"boundary_ranks": [r1, r2], "H1_dimension": 15 - r1 - r2, "H2_dimension": 10 - r2}
    snf = smith_normal_form(boundary[2], domain=sp.ZZ)
    output["nonzero_smith_diagonal_of_boundary2"] = [abs(int(snf[i, i])) for i in range(min(snf.shape)) if snf[i, i]]
    output["Hochster_full_vertex_degree_differences_F2_minus_Q_for_S_mod_I"] = {
        "beta_3_6": output["F2"]["H2_dimension"] - output["Q"]["H2_dimension"],
        "beta_4_6": output["F2"]["H1_dimension"] - output["Q"]["H1_dimension"],
    }
    return output


def continued_fraction_period(d):
    a0 = math.isqrt(d)
    m, denominator, a = 0, 1, a0
    period = []
    while True:
        m = denominator * a - m
        denominator = (d - m * m) // denominator
        a = (a0 + m) // denominator
        period.append(a)
        if a == 2 * a0:
            return period


def pell_examples():
    groups = {}
    for p, q in itertools.combinations(list(sp.primerange(5, 160)), 2):
        if p % 4 != 1 or q % 4 != 1:
            continue
        symbol = 1 if pow(p, (q - 1) // 2, q) == 1 else -1
        period = continued_fraction_period(p * q)
        key = (p % 8, q % 8, symbol)
        record = {"d": p * q, "primes": [p, q], "period": period, "negative_pell_solvable": len(period) % 2 == 1}
        group = groups.setdefault(key, {})
        group.setdefault(record["negative_pell_solvable"], record)
        if len(group) == 2:
            return {"same_signature": {"p_mod8": key[0], "q_mod8": key[1], "legendre_p_over_q": key[2]}, "examples": [group[True], group[False]]}
    raise RuntimeError("No contrasting pair found in the specified finite range")


def main():
    examples = {
        "Hilbert_tangent_seed_1": monomial_hom_dimension([(2, 0, 0), (1, 1, 0), (1, 0, 1), (0, 2, 0), (0, 1, 1), (0, 0, 3)]),
        "Hilbert_tangent_seed_2": monomial_hom_dimension([(2, 0, 0), (1, 1, 0), (1, 0, 1), (0, 2, 0), (0, 0, 2)]),
        "projective_plane_torsion_seed": projective_plane_example(),
        "negative_Pell_seed": pell_examples(),
        "scope": "Known examples and exact finite computations; no new theorem, novelty claim, or general classification.",
    }
    destination = Path(__file__).with_name("research-seed-results.json")
    destination.write_text(json.dumps(examples, ensure_ascii=False, indent=2), encoding="utf-8")
    print(json.dumps(examples, ensure_ascii=False, indent=2))


if __name__ == "__main__":
    main()