"""Independent exact verifier for the level-389 mod-11 Kurihara certificate.

Does not import the producer or NumPy.  Uses scalar sparse elimination,
Merel matrices (with an exact C_q boundary check), balanced Farey parents,
and discrete-log projection to a subgroup of order 11.
"""
import argparse
import csv
import hashlib
import json
from collections import Counter, defaultdict
from fractions import Fraction
from math import gcd
from pathlib import Path
from time import perf_counter

N, P, SIZE = 389, 11, 390
HECKE_PRIMES = (2, 3, 5, 7, 13, 17, 19)
EXPECTED_TRACES = {2: -2, 3: -2, 5: -3, 7: -5, 11: -4,
                   13: -3, 17: -6, 19: 5, 397: 24, 991: -53}
INVERSES = [0] + [pow(x, -1, N) for x in range(1, N)]


class CertificateError(ValueError):
    def __init__(self, step, detail):
        self.step, self.detail = step, detail
        super().__init__(f"{step}: {detail}")


def require(condition, step, detail):
    if not condition:
        raise CertificateError(step, detail)


def projective(c, d):
    c, d = c % N, d % N
    require(c != 0 or d != 0, "projective", "zero bottom row")
    return d * INVERSES[c] % N if c else N


def pair(i):
    return (1, i) if i < N else (0, 1)


def sparse(entries):
    result = {}
    for j, value in entries:
        new = (result.get(j, 0) + value) % P
        if new:
            result[j] = new
        else:
            result.pop(j, None)
    return result


def relation_rows():
    two, three = [], []
    for i in range(SIZE):
        c, d = pair(i)
        two.append(sparse(((i, 1), (projective(d, -c), 1))))
        three.append(sparse(((i, 1), (projective(d, -c-d), 1),
                             (projective(-c-d, c), 1))))
    return two + three


def plus_rows():
    return [sparse(((i, 1), (projective(-pair(i)[0], pair(i)[1]), -1)))
            for i in range(SIZE)]


def matrix(value, rows, columns, step):
    require(isinstance(value, list) and len(value) == rows, step, "wrong row count")
    for i, row in enumerate(value):
        require(isinstance(row, list) and len(row) == columns,
                step, f"wrong column count in row {i}")
        require(all(type(x) is int and 0 <= x < P for x in row),
                step, f"noncanonical field entry in row {i}")
    return value


def sparse_certificate_rows(value, expected, step):
    require(isinstance(value, list) and len(value) == len(expected), step, "wrong row count")
    for i, (entries, wanted) in enumerate(zip(value, expected)):
        require(isinstance(entries, list), step, f"row {i} is not a list")
        indices = []
        for entry in entries:
            require(isinstance(entry, list) and len(entry) == 2,
                    step, f"invalid sparse entry in row {i}")
            j, x = entry
            require(type(j) is int and 0 <= j < SIZE and type(x) is int and 0 < x < P,
                    step, f"invalid coefficient in row {i}")
            indices.append(j)
        require(len(set(indices)) == len(indices), step, f"duplicate index in row {i}")
        require(dict(entries) == wanted, step, f"row {i} differs from defining relation")


def as_sparse(row):
    return {j: x for j, x in enumerate(row) if x}


def sparse_rref(rows, columns):
    """Incremental sparse scalar elimination, then backward substitution."""
    basis = {}
    for source in rows:
        row = dict(source) if isinstance(source, dict) else as_sparse(source)
        while row:
            pivot = min(row)
            if pivot not in basis:
                inverse = pow(row[pivot], -1, P)
                basis[pivot] = {j: x * inverse % P for j, x in row.items()}
                break
            coefficient = row[pivot]
            for j, x in basis[pivot].items():
                new = (row.get(j, 0) - coefficient*x) % P
                if new:
                    row[j] = new
                else:
                    row.pop(j, None)
    pivots = sorted(basis)
    for pivot in reversed(pivots):
        for earlier in pivots:
            if earlier >= pivot:
                break
            coefficient = basis[earlier].get(pivot, 0)
            if coefficient:
                for j, x in basis[pivot].items():
                    new = (basis[earlier].get(j, 0) - coefficient*x) % P
                    if new:
                        basis[earlier][j] = new
                    else:
                        basis[earlier].pop(j, None)
    rref = [[basis[pivot].get(j, 0) for j in range(columns)] for pivot in pivots]
    free = [j for j in range(columns) if j not in basis]
    kernel = [[0]*len(free) for _ in range(columns)]
    for k, j in enumerate(free):
        kernel[j][k] = 1
        for i, pivot in enumerate(pivots):
            kernel[pivot][k] = -rref[i][j] % P
    return rref, pivots, kernel


def multiply(left, right):
    columns = len(right[0]) if right else 0
    output = []
    for source in left:
        source = source if isinstance(source, dict) else as_sparse(source)
        result = [0]*columns
        for j, coefficient in source.items():
            for k, x in enumerate(right[j]):
                result[k] += coefficient*x
        output.append([x % P for x in result])
    return output


def residual(rows, vector, step):
    for i, row in enumerate(rows):
        value = sum(x*vector[j] for j, x in row.items()) % P
        require(value == 0, step, f"nonzero residual {value} in row {i}")


def check_matrix(saved, expected, step):
    matrix(saved, len(expected), len(expected[0]) if expected else 0, step)
    for i, (got, wanted) in enumerate(zip(saved, expected)):
        require(got == wanted, step, f"row {i} differs from independent recomputation")


def stage(saved, prefix, rows, columns):
    rr, pivots, kernel = sparse_rref(rows, columns)
    require(saved[prefix + "pivots"] == pivots, prefix + "pivots", "pivot columns differ")
    check_matrix(saved[prefix + "rref_nonzero"], rr, prefix + "rref")
    return rr, pivots, kernel


def count_curve(ell):
    # Enumerate the left-hand quadratic rather than the producer's discriminants.
    multiplicity = Counter((y*y+y) % ell for y in range(ell))
    return 1 + sum(multiplicity[(x*x*x+x*x-2*x) % ell] for x in range(ell))


def merel_matrices(q):
    """Positive Merel matrices; a+d<=q+1 follows from det>=a+d-1."""
    matrices = []
    for a in range(1, q+1):
        for d in range(1, q+2-a):
            for b in range(a):
                for c in range(d):
                    if a*d-b*c == q:
                        matrices.append((a, b, c, d))
    boundaries = defaultdict(Counter)
    for a, b, c, d in matrices:
        u, v = (a % q, c % q)
        if u == 0 and v == 0:
            u, v = b % q, d % q
        label = "infinity" if v == 0 else u*pow(v, -1, q) % q
        end = "infinity" if c == 0 else Fraction(a, c)
        start = Fraction(b, d)
        boundaries[label][end] += 1
        boundaries[label][start] -= 1
    require(len(boundaries) == q+1, f"merel:{q}", "wrong number of right cosets")
    wanted = {"infinity": 1, Fraction(0): -1}
    for label, boundary in boundaries.items():
        require({cusp: x for cusp, x in boundary.items() if x} == wanted,
                f"merel:{q}", f"condition C_q failed in coset {label}")
    return matrices


def merel_operator(q):
    matrices = merel_matrices(q)
    rows = []
    for i in range(SIZE):
        u, v = pair(i)
        rows.append(sparse((projective(u*a+v*c, u*b+v*d), 1)
                           for a, b, c, d in matrices))
    return rows, len(matrices)


def verify_eigenline(certificate, expected_normalization_scalar=1):
    """Default protocol fixes first nonzero=1; explicit scale is test-only opt-in."""
    require(certificate.get("schema") == "kurihara-finite-eigenline-v1", "schema", "eigenline schema")
    wanted_parameters = {"level": N, "prime": P, "a_invariants": [0, 1, 1, -2, 0],
                         "summation_primes": [397, 991], "primitive_roots": [5, 6]}
    require(certificate["parameters"] == wanted_parameters, "parameters", "unexpected arithmetic inputs")
    require(certificate["generator_order"] == [list(pair(i)) for i in range(SIZE)],
            "generator_order", "generator order differs")
    conventions = certificate["conventions"]
    require(conventions["S"] == [[0, -1], [1, 0]] and conventions["R"] == [[0, -1], [1, -1]],
            "conventions", "S or R matrix changed")
    R, plus = relation_rows(), plus_rows()
    sparse_certificate_rows(certificate["relation_rows"], R, "relation_rows")
    sparse_certificate_rows(certificate["plus_rows"], plus, "plus_rows")
    data = certificate["linear_algebra"]
    vector = data["eigenvector"]
    matrix([vector], 1, SIZE, "eigenvector")
    require(any(vector), "eigenvector", "zero vector")
    # These residual checks deliberately precede artifact hashes/normalization.
    residual(R, vector, "eigenvector:relation_residual")
    residual(plus, vector, "eigenvector:plus_residual")
    rr, pivots, K = stage(data, "relation_", R, SIZE)
    require(data["relation_rank"] == len(pivots) == 325, "relation_rank", "rank must be 325")
    check_matrix(data["manin_kernel"], K, "manin_kernel")
    traces = {}
    for ell, wanted in EXPECTED_TRACES.items():
        count = count_curve(ell)
        trace = ell+1-count
        require(trace == wanted, f"point_count:{ell}", "independent trace differs from expected curve")
        require(certificate["point_counts"][str(ell)] == {"a_ell": trace, "point_count": count},
                f"point_count:{ell}", "saved point count differs")
        traces[ell] = trace
    require(data["hecke_primes"] == [2, 3, 5], "hecke_primes", "isolation primes changed")
    require(set(certificate["hecke_matrices"]) == {str(q) for q in HECKE_PRIMES},
            "hecke_matrices", "missing or additional operator")
    constraints, merel_counts, operators_on_K = [], {}, {}
    operators = {}
    for q in HECKE_PRIMES:
        independent, count = merel_operator(q)
        operators[q] = independent
        merel_counts[str(q)] = count
        saved = matrix(certificate["hecke_matrices"][str(q)], SIZE, SIZE, f"hecke:{q}:matrix")
        on_K = multiply(independent, K)
        saved_on_K = multiply(saved, K)
        require(on_K == saved_on_K, f"hecke:{q}:whole_quotient", "operator differs on 65-dimensional quotient")
        # Confirm the independent operator descends through every Manin relation.
        require(all(not any(row) for row in multiply(R, on_K)),
                f"hecke:{q}:descent", "relation has nonzero image on quotient")
        for i, row in enumerate(independent):
            value = (sum(x*vector[j] for j, x in row.items())-traces[q]*vector[i]) % P
            require(value == 0, f"eigenvector:hecke_residual:{q}", f"row {i} residual {value}")
        operators_on_K[q] = on_K
        if q in (2, 3, 5):
            constraints.extend([[(x-traces[q]*y) % P for x, y in zip(row, K[i])]
                                for i, row in enumerate(on_K)])
    require(multiply(operators[2], operators_on_K[3]) == multiply(operators[3], operators_on_K[2]),
            "hecke:commutator", "T2 and T3 fail to commute on the whole quotient")
    check_matrix(data["hecke_constraints_on_kernel"], constraints, "hecke_constraints_on_kernel")
    er, ep, L = stage(data, "hecke_constraint_", constraints, len(K[0]))
    check_matrix(data["hecke_kernel_coordinates"], L, "hecke_kernel_coordinates")
    B = multiply(K, L)
    check_matrix(data["hecke_eigenbasis"], B, "hecke_eigenbasis")
    require(len(B[0]) == data["hecke_eigenspace_dimension"] == 2,
            "hecke_eigenspace_dimension", "joint eigenspace must have dimension 2")
    plus_constraints = multiply(plus, B)
    check_matrix(data["plus_constraints_on_eigenbasis"], plus_constraints, "plus_constraints_on_eigenbasis")
    pr, pp, C = stage(data, "plus_constraint_", plus_constraints, len(B[0]))
    check_matrix(data["plus_kernel_coordinates"], C, "plus_kernel_coordinates")
    line = multiply(B, C)
    require(len(line[0]) == data["plus_eigenspace_dimension"] == 1,
            "plus_eigenspace_dimension", "geometric plus eigenspace must have dimension 1")
    first = next(i for i, row in enumerate(line) if row[0])
    require(type(expected_normalization_scalar) is int and 1 <= expected_normalization_scalar < P,
            "normalization", "invalid explicit normalization scalar")
    wanted_rule = f"first nonzero coordinate in the fixed generator order equals {expected_normalization_scalar}"
    require(data["normalization_rule"] == wanted_rule, "normalization:rule", "normalization rule changed")
    require(data["normalization_index"] == first, "normalization:index", "first nonzero index differs")
    require(vector[first] == expected_normalization_scalar, "normalization:scalar", "first nonzero entry violates rule")
    scale = expected_normalization_scalar*pow(line[first][0], -1, P) % P
    require(vector == [row[0]*scale % P for row in line], "normalization:vector", "vector differs from normalized line")
    # Every original generator path is also recovered using Farey parents.
    for i in range(SIZE):
        start, end = ((i-1, i), (1, 1)) if i < N else ((1, 1), (1, 0))
        words = sparse([(j, 1) for j in farey_word(*end)] + [(j, -1) for j in farey_word(*start)])
        require(multiply([words], K)[0] == K[i], "farey:generator_path", f"generator {i} failed")
    return {"vector": vector, "kernel": K, "summary": {
        "relation_rank": len(pivots), "manin_dimension": len(K[0]),
        "hecke_constraint_rank": len(ep), "hecke_eigenspace_dimension": len(B[0]),
        "plus_constraint_rank": len(pp), "plus_eigenspace_dimension": len(line[0]),
        "normalization_index": first, "normalization_value": expected_normalization_scalar,
        "merel_matrix_counts": merel_counts,
        "merel_right_coset_boundary_checks": {str(q): q+1 for q in HECKE_PRIMES},
        "whole_quotient_hecke_comparisons": len(HECKE_PRIMES)*SIZE*len(K[0]),
        "hecke_eigenvector_residuals_checked": len(HECKE_PRIMES)*SIZE,
        "farey_generator_paths_checked_on_whole_quotient": SIZE,
        "independent_point_counts": {str(q): q+1-a for q, a in traces.items()},
    }}


def farey_word(a, b):
    """Balanced inverse/Farey-parent path {infinity,a/b}, no continued fractions."""
    if b == 0:
        return ()
    if b < 0:
        a, b = -a, -b
    common = gcd(a, b)
    a, b = a//common, b//common
    a %= b
    word = [projective(1, 0)]
    while b > 1:
        inverse = pow(a, -1, b)
        sign, denominator = (1, inverse) if 2*inverse <= b else (-1, b-inverse)
        numerator = (a*denominator-sign)//b
        require(a*denominator-b*numerator == sign, "farey:determinant", "edge is not unimodular")
        word.append(projective(sign*b, denominator))
        a, b = numerator, denominator
        a %= b
    return tuple(word)


def farey_value(a, b, vector):
    return sum(vector[j] for j in farey_word(a, b)) % P


def prime_factors(n):
    answer = []
    q = 2
    while q*q <= n:
        if n % q == 0:
            answer.append(q)
            while n % q == 0:
                n //= q
        q += 1
    if n > 1:
        answer.append(n)
    return answer


def projected_logs(ell, root):
    require(pow(root, ell-1, ell) == 1 and
            all(pow(root, (ell-1)//q, ell) != 1 for q in prime_factors(ell-1)),
            f"logs:{ell}:root", "primitive-root condition failed")
    exponent = (ell-1)//P
    basis = pow(root, exponent, ell)
    lookup = {pow(basis, j, ell): j for j in range(P)}
    require(len(lookup) == P, f"logs:{ell}:projection", "subgroup is not order 11")
    return [-1] + [lookup[pow(a, exponent, ell)] for a in range(1, ell)]


def verify_logs(certificate):
    require(certificate["modulus"] == P, "logs:modulus", "wrong coefficient field")
    require(set(certificate["tables"]) == {"397", "991"}, "logs:tables", "wrong primes")
    result = {}
    for ell, root in ((397, 5), (991, 6)):
        table = certificate["tables"][str(ell)]
        require(table["root"] == root, f"logs:{ell}:root", "specified primitive root changed")
        values = table["values"]
        wanted = projected_logs(ell, root)
        require(isinstance(values, list) and len(values) == ell,
                f"logs:{ell}:length", "wrong table length")
        for a, (got, expected) in enumerate(zip(values, wanted)):
            require(type(got) is int and got == expected,
                    f"logs:{ell}:entry", f"residue {a} has log {got}, expected {expected}")
        result[ell] = wanted
    return result


def sha256(path):
    h = hashlib.sha256()
    with open(path, "rb") as stream:
        for chunk in iter(lambda: stream.read(1024*1024), b""):
            h.update(chunk)
    return h.hexdigest()


def verify_sum_metadata(metadata, blocks, raw_sum, unit_count):
    require(metadata["n"] == 397*991 and metadata["modulus"] == P,
            "sum:parameters", "wrong conductor or coefficient field")
    require(metadata["unit_count"] == unit_count == 396*990,
            "sum:unit_count", "incorrect unit count")
    require(metadata["raw_product_sum"] == raw_sum,
            "sum:raw_product_sum", "integer sum differs")
    require(metadata["residue"] == raw_sum % P,
            "sum:residue", "mod-11 sum differs")
    require(metadata["blocks"] == blocks, "sum:blocks", "block trace differs")


def verify_sum(results, metadata, vector, logs):
    results = Path(results)
    require(metadata.get("schema") == "kurihara-finite-sum-v1", "sum:schema", "unknown sum schema")
    for key, filename in (("eigenline_sha256", "eigenline_certificate.json"),
                          ("log_tables_sha256", "log_tables.json"),
                          ("terms_csv_sha256", "kurihara_terms.csv")):
        require(metadata[key] == sha256(results/filename), f"sum:hash:{key}", "artifact hash differs")
    n = 397*991
    blocks = [{"first_a": a, "last_a": min(n-1, a+9999), "unit_count": 0,
               "raw_product_sum": 0, "residue": 0} for a in range(1, n, 10000)]
    total, count, reduced_total, word_edges = 0, 0, 0, 0
    with open(results/"kurihara_terms.csv", newline="", encoding="utf-8") as stream:
        reader = csv.reader(stream)
        require(next(reader, None) == ["a", "symbol", "log_397", "log_991", "term"],
                "terms:header", "CSV header differs")
        for a in range(1, n):
            if a % 397 == 0 or a % 991 == 0:
                continue
            row = next(reader, None)
            require(row is not None and len(row) == 5, "terms:row", f"missing or malformed row for a={a}")
            try:
                supplied = [int(x) for x in row]
            except ValueError as exc:
                raise CertificateError("terms:integer", f"noninteger row for a={a}") from exc
            word = farey_word(a, n)
            word_edges += len(word)
            symbol = sum(vector[j] for j in word) % P
            l1, l2 = logs[397][a % 397], logs[991][a % 991]
            raw = symbol*l1*l2
            wanted = [a, symbol, l1, l2, raw % P]
            for name, got, expected in zip(("a", "symbol", "log_397", "log_991", "term"), supplied, wanted):
                require(got == expected, f"terms:{name}", f"a={a}: got {got}, expected {expected}")
            count += 1
            total += raw
            reduced_total += raw % P
            block = blocks[(a-1)//10000]
            block["unit_count"] += 1
            block["raw_product_sum"] += raw
        require(next(reader, None) is None, "terms:extra_row", "unexpected CSV row after final unit")
    for block in blocks:
        block["residue"] = block["raw_product_sum"] % P
    verify_sum_metadata(metadata, blocks, total, count)
    return {"n": n, "modulus": P, "unit_count": count, "raw_product_sum": total,
            "sum_of_reduced_terms": reduced_total, "residue": total % P,
            "block_count": len(blocks), "farey_word_edges_checked": word_edges,
            "csv_fields_checked": 5*count, "all_csv_rows_independently_recomputed": True,
            "archival_literal_6_status": "NOT_REPRODUCED_NORMALIZATION_UNSPECIFIED"}


def verify_all(results):
    started = perf_counter()
    results = Path(results)
    with open(results/"eigenline_certificate.json", encoding="utf-8") as stream:
        eigenline = json.load(stream)
    with open(results/"log_tables.json", encoding="utf-8") as stream:
        log_certificate = json.load(stream)
    with open(results/"kurihara_certificate.json", encoding="utf-8") as stream:
        sum_certificate = json.load(stream)
    eigen = verify_eigenline(eigenline)
    logs = verify_logs(log_certificate)
    summation = verify_sum(results, sum_certificate, eigen["vector"], logs)
    return {"schema": "kurihara-independent-verification-v1", "verdict": "PASS",
            "arithmetic": "Python integers modulo 11; no NumPy or producer imports",
            "independent_methods": {"linear_algebra": "incremental sparse scalar elimination",
                "hecke": "positive Merel matrices plus exact rational C_q boundary checks",
                "modular_paths": "balanced Farey parents using modular inverses",
                "logarithms": "projection to the multiplicative subgroup of order 11",
                "point_counts": "multiplicity table of y^2+y"},
            "linear_algebra": eigen["summary"], "summation": summation,
            "elapsed_seconds": round(perf_counter()-started, 3)}


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--results", type=Path, default=Path(__file__).resolve().parents[1]/"results")
    parser.add_argument("--report", type=Path)
    arguments = parser.parse_args()
    report_path = arguments.report or arguments.results/"independent_verification.json"
    try:
        report = verify_all(arguments.results)
    except (CertificateError, KeyError, ValueError, OSError, TypeError) as exc:
        report = {"schema": "kurihara-independent-verification-v1", "verdict": "FAIL",
                  "failed_step": getattr(exc, "step", "input"), "detail": str(exc)}
    report_path.parent.mkdir(parents=True, exist_ok=True)
    report_path.write_text(json.dumps(report, indent=2, sort_keys=True)+"\n", encoding="utf-8")
    print(json.dumps(report, indent=2, sort_keys=True))
    return 0 if report["verdict"] == "PASS" else 1


if __name__ == "__main__":
    raise SystemExit(main())