#!/usr/bin/env python3
"""Independent source, algebra, and fixture audit for canon claim C-0016."""

from __future__ import annotations

from dataclasses import dataclass
from fractions import Fraction
from hashlib import sha256
from itertools import combinations
import json
from pathlib import Path
import sys

import flint
from flint import arb, arb_mat, ctx
import sympy as sp


if not __debug__:
    raise RuntimeError("optimized Python disables the audit assertions")

assert sys.version_info[:2] == (3, 14)
assert sp.__version__ == "1.14.0"
assert flint.__version__ == "0.9.0"
assert flint.__FLINT_VERSION__ == "3.6.0"


ROOT = Path(__file__).resolve().parent
SOURCES = ROOT / "sources"

SOURCE_HASHES = {
    "milne-abelian-varieties-v2.pdf": (
        "f5ca4e63e5092a4b102daad1470e4cbed"
        "5fe8f82115e3a28c8881e3f67f6aaef"
    ),
    "milne-abelian-varieties-v2.txt": (
        "fb9fcd9960d6df085998d9a79c6784ff"
        "40af49e4edd9dce0c887722785988758"
    ),
    "milne-complex-multiplication.pdf": (
        "0ac337a885d72f6938f15724b897f568f"
        "e350d0d754d0b52370b8f3a01d2e2b5"
    ),
    "milne-complex-multiplication.txt": (
        "c78451f467db58e7c3dd964658a195bc"
        "5837a349ce4ccf3d6aa4b621f65aa165"
    ),
    "milne-jacobian-varieties.pdf": (
        "36c3f09c7462dbbd4ae1f8b81a02bd9f"
        "f84f03c5a346351d7d5d78fc3f173486"
    ),
    "milne-jacobian-varieties.txt": (
        "4f796518221e5faf5c0ddfbf0115c8ca"
        "7d182a1f93ec72e72abb79ae84eee10d"
    ),
}


def file_hash(path: Path) -> str:
    return sha256(path.read_bytes()).hexdigest()


def normalized_text(path: Path) -> str:
    return " ".join(path.read_text(errors="strict").replace("\f", " ").split())


def require_ordered_anchors(text: str, anchors: list[str]) -> None:
    position = 0
    for anchor in anchors:
        found = text.find(anchor, position)
        assert found >= 0, f"missing source anchor: {anchor!r}"
        position = found + len(anchor)


def verify_sources() -> dict[str, object]:
    observed = {}
    for filename, expected in SOURCE_HASHES.items():
        path = SOURCES / filename
        digest = file_hash(path)
        assert digest == expected
        observed[filename] = digest

    av = normalized_text(SOURCES / "milne-abelian-varieties-v2.txt")
    cm = normalized_text(SOURCES / "milne-complex-multiplication.txt")
    jv = normalized_text(SOURCES / "milne-jacobian-varieties.txt")

    require_ordered_anchors(
        av,
        [
            "The next result is very important.",
            "T HEOREM 14.3. The bilinear form",
            "is positive definite",
            "Chapter II Abelian Varieties: Arithmetic",
            "For a variety V over Fq , the Frobenius map",
            "L EMMA 1.2. Let",
            "Rosati involution on End.A/ ˝ Q defined by a polarization",
            "L EMMA 1.3. Let A be an abelian variety over a field k "
            "(not necessarily finite).",
            "We have actually proved the following:",
        ],
    )
    require_ordered_anchors(
        cm,
        [
            "Positive involutions and CM-algebras",
            "P ROPOSITION 1.36 Every finite-dimensional Q-algebra admitting "
            "a positive involution is semisimple.",
            "E XAMPLE 1.38",
            "For a totally real number field F , the identity involution is positive.",
            "For a CM-field E, the involution",
            "P ROPOSITION 1.39 Every finite-dimensional commutative Q-algebra "
            "with positive invo- lution is a product of pairs as in (1.38).",
        ],
    )
    require_ordered_anchors(
        jv,
        [
            "is a principal polarization of J , called the canonical polarization.",
            "C OROLLARY 9.6. For all primes",
            "the map of",
            "H 1 .J; Z",
            "H 1 .C; Z",
            "induced by f is an isomorphism.",
            "11 The zeta function of a curve",
            "the characteristic polynomial of the Frobenius endomorphism",
            "P .X / of degree 2g",
            "T HEOREM 11.1.",
            "The number N of points on C",
            "C OROLLARY 11.4. The zeta function Z.C; t / of C is equal to",
            "R EMARK 11.5. As we saw in (9.6)",
        ],
    )
    characteristic_index = jv.index(
        "the characteristic polynomial of the Frobenius endomorphism"
    )
    corollary_index = jv.index(
        "C OROLLARY 11.4. The zeta function Z.C; t / of C is equal to"
    )
    assert "P .X /" in jv[characteristic_index : characteristic_index + 240]
    assert "P .t /" in jv[corollary_index : corollary_index + 160]

    return {
        "hashes": observed,
        "abelian_varieties_anchors": {
            "rosati_positivity": "Theorem I.14.3, printed pages 61-62",
            "frobenius_rosati": "Lemma II.1.2, printed pages 76-77",
            "embedding_modulus": "Lemma II.1.3 and Remark II.1.4, pages 77-78",
        },
        "cm_cross_check": {
            "semisimplicity": "Proposition I.1.36, printed page 20",
            "commutative_classification": "Proposition I.1.39, printed pages 20-21",
        },
        "jacobian_varieties_2021": {
            "jacobian_polarization": "Summary 6.11, printed pages 22-23",
            "cohomology_identification": "Corollary 9.6, printed page 31",
            "point_count_formula": "Theorem 11.1, printed pages 35-36",
            "printed_convention_mismatch": (
                "The source defines monic P(X)=prod(X-alpha) but Corollary "
                "11.4 prints P(t) as the numerator; the corollary is not "
                "used for that polynomial identity."
            ),
        },
    }


def verify_power_sum_zeta_derivation() -> dict[str, str]:
    t, q = sp.symbols("t q")
    alphas = sp.symbols("alpha_0:4")
    numerator = sp.prod(1 - alpha * t for alpha in alphas)
    zeta_rational = numerator / ((1 - t) * (1 - q * t))
    logarithmic_derivative = sp.factor(
        t * sp.diff(zeta_rational, t) / zeta_rational
    )
    point_count_generating_function = (
        t / (1 - t)
        + q * t / (1 - q * t)
        - sum(
            alpha * t / (1 - alpha * t)
            for alpha in alphas
        )
    )
    assert sp.factor(
        logarithmic_derivative - point_count_generating_function
    ) == 0

    expansion = sp.series(
        point_count_generating_function,
        t,
        0,
        9,
    ).removeO().expand()
    for degree in range(1, 9):
        expected = (
            1 + q**degree
            - sum(alpha**degree for alpha in alphas)
        )
        assert sp.expand(expansion.coeff(t, degree) - expected) == 0

    return {
        "source_input": (
            "Theorem 11.1 applied to C/F_(q^n): "
            "N_n=1+q^n-sum_i alpha_i^n"
        ),
        "exact_map": (
            "t*d/dt log Z_C(t)=t/(1-t)+q*t/(1-q*t)"
            "-sum_i alpha_i*t/(1-alpha_i*t)"
        ),
        "numerator": "P_C(t)=prod_i(1-alpha_i*t)",
        "source_caveat": (
            "Corollary 11.4's printed P(t) is not used because the article "
            "previously defines P(X)=prod_i(X-alpha_i)."
        ),
    }


def verify_symbolic_kernel_bookkeeping() -> dict[str, str]:
    z, ell, beta = sp.symbols("z ell beta", nonzero=True)
    pair = sp.exp(ell * z) * (
        1 - beta * sp.exp(-ell * z)
    ) * (
        1 - sp.exp(-ell * z) / beta
    )
    expected_pair = (
        2 * sp.cosh(ell * z) - beta - 1 / beta
    )
    assert sp.simplify(sp.expand(pair) - expected_pair) == 0

    w = sp.symbols("w")
    assert sp.simplify(
        2 * (sp.cosh(w) - 1) - 4 * sp.sinh(w / 2) ** 2
    ) == 0
    assert sp.simplify(
        2 * (sp.cosh(w) + 1) - 4 * sp.cosh(w / 2) ** 2
    ) == 0

    k = sp.symbols("k", integer=True)
    theta = sp.symbols("theta", real=True)
    generic_zero = sp.I * (2 * sp.pi * k + theta)
    assert sp.trigsimp(
        sp.cosh(generic_zero) - sp.cos(theta)
    ) == 0
    assert sp.trigsimp(sp.sinh(generic_zero)) == sp.I * sp.sin(theta)

    zero_angle = 2 * sp.pi * sp.I * k
    assert sp.simplify(sp.cosh(zero_angle) - 1) == 0
    assert sp.simplify(sp.sinh(zero_angle)) == 0
    assert sp.simplify(sp.cosh(zero_angle)) == 1

    pi_angle = (2 * k + 1) * sp.pi * sp.I
    assert sp.simplify(sp.cosh(pi_angle) + 1) == 0
    assert sp.simplify(sp.sinh(pi_angle)) == 0
    assert sp.simplify(sp.cosh(pi_angle)) == -1

    x, y, gamma = sp.symbols("x y gamma", positive=True)
    atom_x = 1 / (x + gamma**2)
    atom_y = 1 / (y + gamma**2)
    divided_atom = sp.factor(-(atom_x - atom_y) / (x - y))
    assert divided_atom == 1 / (
        (x + gamma**2) * (y + gamma**2)
    )
    zero_atom = sp.factor(-(1 / x - 1 / y) / (x - y))
    assert zero_atom == 1 / (x * y)
    assert sp.simplify(
        -sp.diff(atom_x, x) - 1 / (x + gamma**2) ** 2
    ) == 0
    assert sp.simplify(-sp.diff(1 / x, x) - 1 / x**2) == 0

    lam, mu = sp.symbols("lam mu", positive=True)
    h = sp.Function("h")
    log_lam = lam * h(lam**2)
    log_mu = mu * h(mu**2)
    q_kernel = 4 * (
        mu * log_lam - lam * log_mu
    ) / (mu**2 - lam**2)
    gram_congruence = (
        4
        * lam
        * mu
        * (-(h(lam**2) - h(mu**2)) / (lam**2 - mu**2))
    )
    assert sp.simplify(q_kernel - gram_congruence) == 0

    return {
        "pair_factor": (
            "exp(Lz)(1-beta exp(-Lz))(1-beta^-1 exp(-Lz))="
            "2 cosh(Lz)-beta-beta^-1"
        ),
        "generic_positive_frequencies": (
            "(2*pi*k+theta)/L and (2*pi*(k+1)-theta)/L, k>=0"
        ),
        "theta_zero": (
            "frequency 0 and 2*pi*k/L (k>=1), all doubled per reciprocal pair"
        ),
        "theta_pi": (
            "(2*k+1)*pi/L (k>=0), all doubled per reciprocal pair"
        ),
        "canonical_product": (
            "X(z)=C*z^m0*product_(gamma>0)"
            "(1+z^2/gamma^2)^m_gamma"
        ),
        "log_derivative": (
            "H(x)=m0/x+2*sum_(gamma>0)m_gamma/(x+gamma^2)"
        ),
        "confluent_value": (
            "M(x,x)=-H'(x)=m0/x^2+2*sum m_gamma/(x+gamma^2)^2"
        ),
        "congruence": "Q(lambda_i,lambda_j)=(2 lambda_i) M_ij (2 lambda_j)",
    }


@dataclass(frozen=True)
class FiniteField:
    p: int
    modulus: tuple[int, ...]

    def __post_init__(self) -> None:
        assert self.modulus[-1] % self.p == 1
        assert len(self.modulus) >= 2
        variable = sp.symbols("t")
        polynomial = sp.Poly(
            sum(
                (coefficient % self.p) * variable**degree
                for degree, coefficient in enumerate(self.modulus)
            ),
            variable,
            modulus=self.p,
        )
        assert polynomial.degree() == self.degree
        assert polynomial.is_irreducible

    @property
    def degree(self) -> int:
        return len(self.modulus) - 1

    @property
    def order(self) -> int:
        return self.p**self.degree

    def digits(self, value: int) -> list[int]:
        output = []
        for _ in range(self.degree):
            output.append(value % self.p)
            value //= self.p
        return output

    def encode(self, coefficients: tuple[int, ...]) -> int:
        assert len(coefficients) <= self.degree
        return sum(
            (coefficient % self.p) * self.p**index
            for index, coefficient in enumerate(coefficients)
        )

    def add(self, left: int, right: int) -> int:
        return sum(
            ((a + b) % self.p) * self.p**index
            for index, (a, b) in enumerate(
                zip(self.digits(left), self.digits(right), strict=True)
            )
        )

    def neg(self, value: int) -> int:
        return sum(
            ((-coefficient) % self.p) * self.p**index
            for index, coefficient in enumerate(self.digits(value))
        )

    def sub(self, left: int, right: int) -> int:
        return self.add(left, self.neg(right))

    def mul(self, left: int, right: int) -> int:
        left_coefficients = self.digits(left)
        right_coefficients = self.digits(right)
        product = [0] * (2 * self.degree - 1)
        for left_degree, left_coefficient in enumerate(left_coefficients):
            for right_degree, right_coefficient in enumerate(right_coefficients):
                product[left_degree + right_degree] = (
                    product[left_degree + right_degree]
                    + left_coefficient * right_coefficient
                ) % self.p

        for degree in range(
            2 * self.degree - 2,
            self.degree - 1,
            -1,
        ):
            leading = product[degree] % self.p
            for lower_degree in range(self.degree):
                product[degree - self.degree + lower_degree] = (
                    product[degree - self.degree + lower_degree]
                    - leading * self.modulus[lower_degree]
                ) % self.p

        return sum(
            product[index] * self.p**index
            for index in range(self.degree)
        )

    def pow(self, value: int, exponent: int) -> int:
        result = self.encode((1,))
        while exponent:
            if exponent & 1:
                result = self.mul(result, value)
            value = self.mul(value, value)
            exponent //= 2
        return result


def evaluate_sparse_polynomial(
    field: FiniteField,
    value: int,
    terms: tuple[tuple[int, int], ...],
) -> int:
    result = 0
    for exponent, coefficient in terms:
        result = field.add(
            result,
            field.mul(
                field.encode((coefficient,)),
                field.pow(value, exponent),
            ),
        )
    return result


def count_odd_hyperelliptic_points(
    field: FiniteField,
    terms: tuple[tuple[int, int], ...],
    multiplier: int | None = None,
) -> int:
    one = field.encode((1,))
    minus_one = field.encode((-1,))
    total = 1
    for value in range(field.order):
        right = evaluate_sparse_polynomial(field, value, terms)
        if multiplier is not None:
            right = field.mul(multiplier, right)
        if right == 0:
            total += 1
            continue
        character = field.pow(right, (field.order - 1) // 2)
        if character == one:
            total += 2
        else:
            assert character == minus_one
    return total


def assert_squarefree(
    p: int,
    terms: tuple[tuple[int, int], ...],
) -> None:
    variable = sp.symbols("x")
    polynomial = sp.Poly(
        sum(coefficient * variable**exponent for exponent, coefficient in terms),
        variable,
        modulus=p,
    )
    assert sp.gcd(polynomial, polynomial.diff()).degree() == 0


def numerator_from_power_sums(power_sums: list[int]) -> list[int]:
    elementary: list[Fraction] = [Fraction(1)]
    for degree in range(1, len(power_sums) + 1):
        value = sum(
            (-1) ** (index - 1)
            * elementary[degree - index]
            * power_sums[index - 1]
            for index in range(1, degree + 1)
        ) / degree
        elementary.append(value)
    assert all(value.denominator == 1 for value in elementary)
    return [
        int(((-1) ** degree) * elementary[degree])
        for degree in range(len(elementary))
    ]


def verify_curve_fixtures() -> dict[str, object]:
    elliptic_generic = ((3, 1), (1, 1), (0, 1))
    elliptic_trace_zero = ((3, 1), (1, -1))
    genus_two = ((5, 1), (1, -1))
    assert_squarefree(5, elliptic_generic)
    assert_squarefree(3, elliptic_trace_zero)
    assert_squarefree(3, genus_two)
    assert_squarefree(5, genus_two)

    f5 = FiniteField(5, (0, 1))
    f25 = FiniteField(5, (2, 0, 1))
    generic_counts = [
        count_odd_hyperelliptic_points(f5, elliptic_generic),
        count_odd_hyperelliptic_points(f25, elliptic_generic),
    ]
    assert generic_counts == [9, 27]
    generic_power_sums = [
        5**degree + 1 - count
        for degree, count in enumerate(generic_counts, start=1)
    ]
    generic_numerator = numerator_from_power_sums(generic_power_sums)
    assert generic_numerator == [1, 3, 5]

    f9 = FiniteField(3, (1, 0, 1))
    base_count = count_odd_hyperelliptic_points(
        f9,
        elliptic_trace_zero,
    )
    nonsquare = f9.encode((1, 1))
    assert f9.pow(nonsquare, 4) == f9.encode((-1,))
    twist_count = count_odd_hyperelliptic_points(
        f9,
        elliptic_trace_zero,
        multiplier=nonsquare,
    )
    assert (base_count, twist_count) == (16, 4)
    base_trace = 9 + 1 - base_count
    twist_trace = 9 + 1 - twist_count
    base_numerator = [1, -base_trace, 9]
    twist_numerator = [1, -twist_trace, 9]
    assert base_numerator == [1, 6, 9]
    assert twist_numerator == [1, -6, 9]

    extension_moduli = {
        1: (0, 1),
        2: (1, 0, 1),
        3: (1, 0, 2, 1),
        4: (2, 1, 0, 0, 1),
    }
    genus_two_counts = [
        count_odd_hyperelliptic_points(
            FiniteField(3, extension_moduli[degree]),
            genus_two,
        )
        for degree in range(1, 5)
    ]
    assert genus_two_counts == [4, 6, 28, 110]
    genus_two_power_sums = [
        3**degree + 1 - count
        for degree, count in enumerate(genus_two_counts, start=1)
    ]
    genus_two_numerator = numerator_from_power_sums(genus_two_power_sums)
    assert genus_two_numerator == [1, 0, -2, 0, 9]

    exceptional_genus_two_moduli = {
        1: (0, 1),
        2: (1, 1, 1),
        3: (1, 0, 1, 1),
        4: (1, 0, 1, 1, 1),
    }
    exceptional_genus_two_counts = [
        count_odd_hyperelliptic_points(
            FiniteField(5, exceptional_genus_two_moduli[degree]),
            genus_two,
        )
        for degree in range(1, 5)
    ]
    assert exceptional_genus_two_counts == [6, 6, 126, 526]
    exceptional_genus_two_power_sums = [
        5**degree + 1 - count
        for degree, count in enumerate(
            exceptional_genus_two_counts,
            start=1,
        )
    ]
    exceptional_genus_two_numerator = numerator_from_power_sums(
        exceptional_genus_two_power_sums
    )
    assert exceptional_genus_two_numerator == [1, 0, -10, 0, 25]

    z, ell = sp.symbols("z ell")
    generic_centered = sp.exp(ell * z) * (
        1
        + 3 / sp.sqrt(5) * sp.exp(-ell * z)
        + sp.exp(-2 * ell * z)
    )
    assert sp.simplify(
        generic_centered
        - (2 * sp.cosh(ell * z) + 3 / sp.sqrt(5))
    ) == 0
    base_centered = sp.exp(ell * z) * (
        1 + sp.exp(-ell * z)
    ) ** 2
    twist_centered = sp.exp(ell * z) * (
        1 - sp.exp(-ell * z)
    ) ** 2
    assert sp.simplify(
        base_centered - 2 * (sp.cosh(ell * z) + 1)
    ) == 0
    assert sp.simplify(
        twist_centered - 2 * (sp.cosh(ell * z) - 1)
    ) == 0
    genus_two_centered = sp.exp(2 * ell * z) * (
        1
        - sp.Rational(2, 3) * sp.exp(-2 * ell * z)
        + sp.exp(-4 * ell * z)
    )
    assert sp.simplify(
        genus_two_centered
        - (2 * sp.cosh(2 * ell * z) - sp.Rational(2, 3))
    ) == 0
    exceptional_genus_two_centered = sp.exp(2 * ell * z) * (
        1
        - 2 * sp.exp(-2 * ell * z)
        + sp.exp(-4 * ell * z)
    )
    assert sp.simplify(
        exceptional_genus_two_centered
        - (4 * sp.sinh(ell * z) ** 2).rewrite(sp.exp)
    ) == 0

    return {
        "generic_elliptic_F5": {
            "curve": "y^2=x^3+x+1",
            "counts_N1_N2": generic_counts,
            "numerator": generic_numerator,
            "centered_X": "2*cosh(log(5)*z)+3/sqrt(5)",
        },
        "supersingular_elliptic_F9": {
            "curve": "y^2=x^3-x (base change from F3)",
            "count_N1": base_count,
            "numerator": base_numerator,
            "centered_X": "2*(cosh(log(9)*z)+1)",
            "exceptional_angle": "theta=pi; every frequency is doubled",
        },
        "supersingular_twist_F9": {
            "curve": "y^2=(1+t)*(x^3-x), t^2+1=0",
            "count_N1": twist_count,
            "numerator": twist_numerator,
            "centered_X": "2*(cosh(log(9)*z)-1)",
            "exceptional_angle": "theta=0; m0=2 and every nonzero frequency is doubled",
        },
        "genus_two_F3": {
            "curve": "y^2=x^5-x",
            "counts_N1_to_N4": genus_two_counts,
            "power_sums_S1_to_S4": genus_two_power_sums,
            "numerator": genus_two_numerator,
            "centered_X": "2*cosh(2*log(3)*z)-2/3",
        },
        "exceptional_genus_two_F5": {
            "curve": "y^2=x^5-x",
            "counts_N1_to_N4": exceptional_genus_two_counts,
            "power_sums_S1_to_S4": exceptional_genus_two_power_sums,
            "numerator": exceptional_genus_two_numerator,
            "centered_X": "4*sinh(log(5)*z)^2",
            "exceptional_angles": (
                "theta=0 and theta=pi, each reciprocal pair doubled; "
                "m0=2 and every positive frequency has multiplicity 2"
            ),
        },
    }


def arb_polynomial_data(
    coefficients: list[int],
    value: arb,
) -> tuple[arb, arb, arb]:
    polynomial = arb(0)
    first = arb(0)
    second = arb(0)
    for degree, coefficient in enumerate(coefficients):
        polynomial += coefficient * value**degree
        if degree >= 1:
            first += degree * coefficient * value ** (degree - 1)
        if degree >= 2:
            second += (
                degree
                * (degree - 1)
                * coefficient
                * value ** (degree - 2)
            )
    return polynomial, first, second


def descended_value_and_confluent(
    q: int,
    genus: int,
    coefficients: list[int],
    x: arb,
) -> tuple[arb, arb]:
    radius = x.sqrt()
    ell = arb(q).log()
    u = (-ell * radius).exp() / arb(q).sqrt()
    polynomial, first, second = arb_polynomial_data(coefficients, u)
    assert not polynomial.contains(0)
    logarithmic = (
        genus * ell - ell * u * first / polynomial
    )
    f_prime = (
        first / polynomial
        + u
        * (second * polynomial - first**2)
        / polynomial**2
    )
    logarithmic_derivative = ell**2 * u * f_prime
    descended = logarithmic / radius
    negative_x_derivative = (
        logarithmic - radius * logarithmic_derivative
    ) / (2 * radius**3)
    return descended, negative_x_derivative


def principal_minors(matrix: list[list[arb]]) -> dict[str, str]:
    size = len(matrix)
    output = {}
    for order in range(1, size + 1):
        for indices in combinations(range(size), order):
            submatrix = arb_mat(
                [
                    [matrix[row][column] for column in indices]
                    for row in indices
                ]
            )
            determinant = submatrix.det()
            assert determinant > 0
            output[",".join(str(index) for index in indices)] = str(determinant)
    return output


def verify_loewner_fixtures() -> dict[str, object]:
    ctx.prec = 384
    lambdas = [arb(1) / 4, arb(3) / 4, arb(5) / 4]
    xs = [value**2 for value in lambdas]
    fixtures = {
        "generic_elliptic_F5": (5, 1, [1, 3, 5]),
        "supersingular_elliptic_F9": (9, 1, [1, 6, 9]),
        "supersingular_twist_F9": (9, 1, [1, -6, 9]),
        "genus_two_F3": (3, 2, [1, 0, -2, 0, 9]),
        "exceptional_genus_two_F5": (5, 2, [1, 0, -10, 0, 25]),
    }
    output = {}
    for name, (q, genus, coefficients) in fixtures.items():
        values = [
            descended_value_and_confluent(
                q,
                genus,
                coefficients,
                x,
            )
            for x in xs
        ]
        matrix = []
        for row, x in enumerate(xs):
            matrix_row = []
            for column, y in enumerate(xs):
                if row == column:
                    entry = values[row][1]
                else:
                    entry = -(
                        values[row][0] - values[column][0]
                    ) / (x - y)
                assert entry > 0
                matrix_row.append(entry)
            matrix.append(matrix_row)

        q_matrix = [
            [
                4 * lambdas[row] * lambdas[column] * matrix[row][column]
                for column in range(len(lambdas))
            ]
            for row in range(len(lambdas))
        ]
        matrix_det = arb_mat(matrix).det()
        q_det = arb_mat(q_matrix).det()
        scale = arb(1)
        for value in lambdas:
            scale *= (2 * value) ** 2
        assert matrix_det > 0
        assert q_det > 0
        assert (q_det - scale * matrix_det).contains(0)

        output[name] = {
            "lambda_nodes": ["1/4", "3/4", "5/4"],
            "negative_loewner_principal_minors": principal_minors(matrix),
            "Q_principal_minors": principal_minors(q_matrix),
            "determinant_congruence_scale": str(scale),
            "negative_loewner_det": str(matrix_det),
            "Q_det": str(q_det),
        }
    return output


def main() -> None:
    results = {
        "sources": verify_sources(),
        "power_sum_zeta_derivation": verify_power_sum_zeta_derivation(),
        "symbolic_kernel_bookkeeping": verify_symbolic_kernel_bookkeeping(),
        "curve_fixtures": verify_curve_fixtures(),
        "loewner_fixtures": verify_loewner_fixtures(),
        "scope": {
            "frobenius_convention": (
                "q-power scheme Frobenius endomorphism; no geometric/arithmetic "
                "Galois-Frobenius label is used"
            ),
            "kernel_domain": (
                "all positive lambda nodes, with repeated nodes interpreted "
                "confluently; this strengthens the prior lambda>1/2 wording"
            ),
            "number_field_consequence": "none",
        },
    }
    output_path = ROOT / "results.json"
    output_path.write_text(json.dumps(results, indent=2) + "\n")
    print("Milne source chain pinned and anchored")
    print("zero multiplicities and Loewner congruence checked symbolically")
    print("five exact finite-field fixtures reconstructed")
    print("all 3x3 fixture Loewner and Q principal minors certified positive")


if __name__ == "__main__":
    main()
