"""Exact finite arithmetic for 389a1; Python standard library only.

No modular-symbol, Selmer, Sha, p-adic-height, or analytic-BSD computation
is performed by this module. Points use the short model stated in the report.
"""
from dataclasses import dataclass
from math import isqrt


def is_prime(n):
    return n >= 2 and all(n % d for d in range(2, isqrt(n) + 1))


@dataclass(frozen=True)
class Curve:
    ell: int
    a: int = -3024
    b: int = 46224

    def __post_init__(self):
        if self.ell <= 3 or not is_prime(self.ell):
            raise ValueError("An odd prime greater than three is required")
        if (4 * self.a**3 + 27 * self.b**2) % self.ell == 0:
            raise ValueError("Singular reduction")

    def on_curve(self, point):
        return point is None or (
            point[1]**2 - point[0]**3 - self.a * point[0] - self.b
        ) % self.ell == 0

    def neg(self, point):
        return None if point is None else (point[0] % self.ell, -point[1] % self.ell)

    def add(self, point, other):
        if point is None:
            return other
        if other is None:
            return point
        q = self.ell
        x, y = point
        u, v = other
        if (x - u) % q == 0 and (y + v) % q == 0:
            return None
        if (x - u) % q == 0:
            slope = (3 * x*x + self.a) * pow(2*y, -1, q) % q
        else:
            slope = (v-y) * pow(u-x, -1, q) % q
        rx = (slope*slope - x-u) % q
        return rx, (slope*(x-rx)-y) % q

    def mul(self, n, point):
        if n < 0:
            return self.mul(-n, self.neg(point))
        result = None
        while n:
            if n & 1:
                result = self.add(result, point)
            point = self.add(point, point)
            n //= 2
        return result

    def points(self):
        # Independent enumeration of both coordinates, using a table of squares.
        squares = {}
        for y in range(self.ell):
            squares.setdefault(y*y % self.ell, []).append(y)
        result = [None]
        for x in range(self.ell):
            for y in squares.get((x**3 + self.a*x + self.b) % self.ell, []):
                result.append((x, y))
        return result

    def count(self):
        # Character-sum count, independent of group operations.
        total = 1
        for x in range(self.ell):
            value = (x**3 + self.a*x + self.b) % self.ell
            symbol = pow(value, (self.ell-1)//2, self.ell)
            total += 1 + (0 if value == 0 else 1 if symbol == 1 else -1)
        return total

    def orbit(self, point):
        result = []
        current = None
        # Hasse upper bound supplies a finite failure condition.
        for _ in range(self.ell + 2 + isqrt(4*self.ell)):
            result.append(current)
            current = self.add(current, point)
            if current is None:
                return result
        raise ArithmeticError("Orbit did not close within Hasse bound")


def rank_mod(matrix, p=11):
    a = [[v % p for v in row] for row in matrix]
    r = 0
    for c in range(len(a[0]) if a else 0):
        pivot = next((i for i in range(r, len(a)) if a[i][c]), None)
        if pivot is None:
            continue
        a[r], a[pivot] = a[pivot], a[r]
        factor = pow(a[r][c], -1, p)
        a[r] = [v*factor % p for v in a[r]]
        for i in range(len(a)):
            if i != r:
                factor = a[i][c]
                a[i] = [(v-factor*w) % p for v, w in zip(a[i], a[r])]
        r += 1
        if r == len(a):
            break
    return r


def certificate():
    p = 11
    # Minimal generalized Weierstrass model; safe also in characteristic two.
    a1, a2, a3, a4, a6 = 0, 1, 1, -2, 0
    b2 = a1*a1 + 4*a2
    b4 = a1*a3 + 2*a4
    b6 = a3*a3 + 4*a6
    b8 = a1*a1*a6 + 4*a2*a6 - a1*a3*a4 + a2*a3*a3 - a4*a4
    discriminant = -b2*b2*b8 - 8*b4**3 - 27*b6*b6 + 9*b2*b4*b6
    c4 = b2*b2 - 24*b4
    counts_minimal = {}
    for q in (2, 5, 11):
        counts_minimal[str(q)] = 1 + sum(
            (y*y+a1*x*y+a3*y-x**3-a2*x*x-a4*x-a6) % q == 0
            for x in range(q) for y in range(q)
        )
    rows = []
    data = []
    for ell in (397, 991):
        curve = Curve(ell)
        P, Q = (12, 108), (48, 108)
        if not curve.on_curve(P) or not curve.on_curve(Q):
            raise ArithmeticError("Input point is off curve")
        points, order = curve.points(), curve.count()
        if len(points) != order:
            raise ArithmeticError("Independent point counts disagree")
        orbit = curve.orbit(P)
        if set(orbit) != set(points):
            raise ArithmeticError("P does not generate the full group")
        m = order // p
        if order % p or m % p == 0:
            raise ArithmeticError("Expected exactly one factor of eleven")
        pP, pQ = curve.mul(m, P), curve.mul(m, Q)
        small_orbit = curve.orbit(pP)
        ratio = small_orbit.index(pQ)
        rows.append([1, ratio])
        data.append({
            "ell": ell, "group_order": order,
            "a_ell": ell+1-order, "order_of_P": len(orbit),
            "Q_discrete_log_base_P": orbit.index(Q),
            "prime_to_11_multiplier": m,
            "multiplied_P": pP, "multiplied_Q": pQ,
            "order_of_multiplied_P": len(small_orbit),
            "Q_norm_coordinate": ratio,
            "independent_counts_agree": True,
            "P_orbit_equals_all_points": True,
        })
    det = (rows[0][0]*rows[1][1]-rows[0][1]*rows[1][0]) % p
    # A lattice basis is independently checked against both congruences.
    basis = [(5742, -22), (-254980, 1045)]
    congruences = []
    for a, b in basis:
        congruences.append([(a+244*b) % 374, (a+356*b) % 1045])
    index = abs(basis[0][0]*basis[1][1]-basis[1][0]*basis[0][1])
    # Enumerate the proposed two-dimensional source, not the full Selmer group.
    images = {
        ((a+2*b) % p, (a+4*b) % p)
        for a in range(p) for b in range(p)
    }
    curve11 = Curve(11)
    return {
        "coefficient_field": "F_11",
        "curve": "y^2+y=x^3+x^2-2x",
        "short_model": "Y^2=X^3-3024X+46224",
        "coordinate_change": {"X": "36*x+12", "Y": "216*y+108"},
        "minimal_invariants": {"b2": b2, "b4": b4, "b6": b6, "b8": b8,
                               "c4": c4, "discriminant": discriminant,
                               "discriminant_is_prime": is_prime(discriminant)},
        "minimal_model_point_counts": counts_minimal,
        "residual_irreducibility_witness": {
            "frobenius_prime": 2,
            "a_2": 3-counts_minimal["2"],
            "characteristic_polynomial_mod_11": [1, 2, 2],
            "discriminant_mod_11": 7,
            "square_residues_mod_11": sorted({x*x % 11 for x in range(11)}),
            "roots_mod_11": [x for x in range(11) if (x*x+2*x+2) % 11 == 0],
            "note": "Transvection and subgroup arguments are proved in the manuscript, not computed here.",
        },
        "good_at_11": {"group_order": curve11.count(), "a_11": 12-curve11.count()},
        "local_results": data, "matrix": rows, "determinant_mod_11": det,
        "source_images_count": len(images),
        "mixed_polynomial_coefficients": {"X^2": 0, "X*Y": det, "Y^2": 0},
        "simultaneous_reduction_lattice": {
            "basis": basis, "congruence_residues": congruences,
            "index": index, "target_order": 374*1045,
        },
        "falsifying_witnesses": {
            "nonzero_mod_121_but_zero_mod_11": {"integer": 11, "mod_121": 11, "mod_11": 0},
            "duplicate_row_determinant": 0,
            "abstract_selmer_countermodel": {
                "matrix": [[1, 2, 0, 0], [1, 4, 0, 0]],
                "source_dimension": 4, "image_dimension": 2, "kernel_dimension": 2,
                "elliptic_curve_realization_claimed": False,
            },
            "same_special_fiber_different_deformation": {
                "fiber_differential": [[0, 0], [0, 0]],
                "lift_0": [["0", "0"], ["0", "0"]],
                "lift_1": [["X", "2*X"], ["Y", "4*Y"]],
                "first_bockstein_zero": True, "second_bockstein_zero": False,
                "determinant_in_associated_graded_degree_2": "2*X*Y",
                "determinant_inside_square_zero_ring_itself": "0",
            },
        },
        "not_computed": [
            "canonical modular symbols or Kurihara number",
            "full Selmer group", "Sha", "global Selmer-complex regulator",
            "complex L-value or height", "formal proof assistant verification",
        ],
    }