verify_discriminant_attack.py

#!/usr/bin/env python3
"""Exact checks for the second Bethe--Calogero structural attack.

The script has no third-party dependencies.  It verifies:

1. the all-orders centered-average factorization of the cubic-divided
   logarithmic theta block through differential order twelve;
2. the displayed invariant coefficients through differential order eight;
3. the rational elliptic-curve formula for the derivative of the Cayley
   root block; and
4. trace discriminants for the split, persistently nonreduced, and mixed
   corank-one model algebras.
"""

from fractions import Fraction
from math import factorial


Poly = dict[tuple[int, int], Fraction]
Series = dict[int, Poly]


def clean(poly: Poly) -> Poly:
    return {monomial: coefficient for monomial, coefficient in poly.items() if coefficient}


def add(left: Poly, right: Poly) -> Poly:
    answer = dict(left)
    for monomial, coefficient in right.items():
        answer[monomial] = answer.get(monomial, Fraction(0)) + coefficient
    return clean(answer)


def scale(poly: Poly, scalar: int | Fraction) -> Poly:
    scalar = Fraction(scalar)
    return clean({monomial: scalar * coefficient for monomial, coefficient in poly.items()})


def multiply(left: Poly, right: Poly) -> Poly:
    answer: Poly = {}
    for (i, j), a in left.items():
        for (k, ell), b in right.items():
            monomial = (i + k, j + ell)
            answer[monomial] = answer.get(monomial, Fraction(0)) + a * b
    return clean(answer)


def power(poly: Poly, exponent: int) -> Poly:
    answer: Poly = {(0, 0): Fraction(1)}
    base = poly
    n = exponent
    while n:
        if n & 1:
            answer = multiply(answer, base)
        base = multiply(base, base)
        n >>= 1
    return answer


def multiply_series(left: Series, right: Series, maximum: int) -> Series:
    answer: Series = {}
    for i, a in left.items():
        for j, b in right.items():
            if i + j <= maximum:
                answer[i + j] = add(answer.get(i + j, {}), multiply(a, b))
    return answer


d1: Poly = {(1, 0): Fraction(1)}
d2: Poly = {(0, 1): Fraction(1)}
d3 = scale(add(d1, d2), -1)
deltas = (d1, d2, d3)
e2 = add(add(multiply(d1, d2), multiply(d2, d3)), multiply(d3, d1))
delta = multiply(multiply(d1, d2), d3)


def power_sum(exponent: int) -> Poly:
    answer: Poly = {}
    for flavor in deltas:
        answer = add(answer, power(flavor, exponent))
    return answer


maximum_order = 12
average_operator: Series = {0: {(0, 0): Fraction(1)}}
for flavor in deltas:
    factor: Series = {}
    for k in range(maximum_order // 2 + 1):
        differential_order = 2 * k
        coefficient = Fraction(1, (2 ** differential_order) * factorial(differential_order + 1))
        factor[differential_order] = scale(power(flavor, differential_order), coefficient)
    average_operator = multiply_series(average_operator, factor, maximum_order)

# The coefficient of D^(n-3) in the product of the three sinhc operators
# must equal 2*p_n/(n!*delta), where p_n=sum_a Delta_a^n.
for differential_order in range(0, maximum_order + 1, 2):
    n = differential_order + 3
    left = multiply(average_operator[differential_order], delta)
    right = scale(power_sum(n), Fraction(2, factorial(n)))
    assert left == right, (differential_order, left, right)

expected_coefficients = {
    0: {(0, 0): Fraction(1)},
    2: scale(e2, Fraction(-1, 12)),
    4: scale(power(e2, 2), Fraction(1, 360)),
    6: add(
        scale(power(e2, 3), Fraction(-1, 20160)),
        scale(power(delta, 2), Fraction(1, 60480)),
    ),
    8: scale(
        add(power(e2, 4), scale(multiply(e2, power(delta, 2)), -1)),
        Fraction(1, 1814400),
    ),
}
for differential_order, expected in expected_coefficients.items():
    assert average_operator[differential_order] == expected, differential_order


def determinant_bareiss(matrix: list[list[int]]) -> int:
    """Fraction-free determinant over the integers."""

    size = len(matrix)
    work = [row[:] for row in matrix]
    sign = 1
    previous = 1
    for pivot_index in range(size - 1):
        if work[pivot_index][pivot_index] == 0:
            swap = next(
                (row for row in range(pivot_index + 1, size) if work[row][pivot_index]),
                None,
            )
            if swap is None:
                return 0
            work[pivot_index], work[swap] = work[swap], work[pivot_index]
            sign *= -1
        pivot = work[pivot_index][pivot_index]
        for i in range(pivot_index + 1, size):
            for j in range(pivot_index + 1, size):
                numerator = work[i][j] * pivot - work[i][pivot_index] * work[pivot_index][j]
                assert numerator % previous == 0
                work[i][j] = numerator // previous
        previous = pivot
    return sign * work[-1][-1]


def split_trace_matrix(mu: int, parameter: int) -> list[list[int]]:
    """Trace matrix in 1,x,...,x^(mu-1) for x^mu=parameter."""

    matrix: list[list[int]] = []
    for i in range(mu):
        row: list[int] = []
        for j in range(mu):
            exponent = i + j
            if exponent == 0:
                row.append(mu)
            elif exponent == mu:
                row.append(mu * parameter)
            else:
                row.append(0)
        matrix.append(row)
    return matrix


for mu in range(2, 10):
    for parameter in (2, 3):
        observed = determinant_bareiss(split_trace_matrix(mu, parameter))
        expected = (
            (-1) ** (mu * (mu - 1) // 2)
            * (mu ** mu)
            * ((-parameter) ** (mu - 1))
        )
        assert observed == expected, (mu, parameter, observed, expected)

# A constant length-mu fat point has trace matrix diag(mu,0,...,0), hence
# zero discriminant.  The mixed flat model x(x-t)^2 has two generic supports
# but is still nonreduced and also has zero trace discriminant.
for mu in range(2, 10):
    assert determinant_bareiss(split_trace_matrix(mu, 0)) == 0

mixed_parameter = 2
mixed_moments = [3] + [2 * mixed_parameter**k for k in range(1, 5)]
mixed_trace_matrix = [
    [mixed_moments[i + j] for j in range(3)]
    for i in range(3)
]
assert determinant_bareiss(mixed_trace_matrix) == 0

# If R=(1-T)/(1+T), T=c*y/(1+a*x), then
# d(log R)/dz = -2*c*N/D on y^2=4*x^3-g2*x-g3.  Check the reduction of N as
# an exact polynomial identity in the independent symbols (x,g2,g3,a).
Poly4 = dict[tuple[int, int, int, int], Fraction]


def add4(left: Poly4, right: Poly4) -> Poly4:
    answer = dict(left)
    for monomial, coefficient in right.items():
        answer[monomial] = answer.get(monomial, Fraction(0)) + coefficient
    return {monomial: coefficient for monomial, coefficient in answer.items() if coefficient}


def scale4(poly: Poly4, scalar: int | Fraction) -> Poly4:
    return {
        monomial: Fraction(scalar) * coefficient
        for monomial, coefficient in poly.items()
        if coefficient
    }


def multiply4(left: Poly4, right: Poly4) -> Poly4:
    answer: Poly4 = {}
    for left_monomial, left_coefficient in left.items():
        for right_monomial, right_coefficient in right.items():
            monomial = tuple(
                left_exponent + right_exponent
                for left_exponent, right_exponent in zip(left_monomial, right_monomial)
            )
            answer[monomial] = (
                answer.get(monomial, Fraction(0))
                + left_coefficient * right_coefficient
            )
    return {monomial: coefficient for monomial, coefficient in answer.items() if coefficient}


one4: Poly4 = {(0, 0, 0, 0): Fraction(1)}
x4: Poly4 = {(1, 0, 0, 0): Fraction(1)}
g2_4: Poly4 = {(0, 1, 0, 0): Fraction(1)}
g3_4: Poly4 = {(0, 0, 1, 0): Fraction(1)}
a4: Poly4 = {(0, 0, 0, 1): Fraction(1)}
x_squared4 = multiply4(x4, x4)
x_cubed4 = multiply4(x_squared4, x4)
y_squared4 = add4(
    add4(scale4(x_cubed4, 4), scale4(multiply4(g2_4, x4), -1)),
    scale4(g3_4, -1),
)
y_prime4 = add4(scale4(x_squared4, 6), scale4(g2_4, Fraction(-1, 2)))
unreduced4 = add4(
    multiply4(y_prime4, add4(one4, multiply4(a4, x4))),
    scale4(multiply4(a4, y_squared4), -1),
)
reduced4 = add4(
    y_prime4,
    multiply4(
        a4,
        add4(
            add4(scale4(x_cubed4, 2), scale4(multiply4(g2_4, x4), Fraction(1, 2))),
            g3_4,
        ),
    ),
)
assert unreduced4 == reduced4

print("exact centered-average operator identity checked through D^12")
print("A_Delta = product_a sinh(Delta_a D/2)/(Delta_a D/2)")
print("coefficients:")
print("  D^0: 1")
print("  D^2: -e2/12")
print("  D^4: e2^2/360")
print("  D^6: -e2^3/20160 + delta^2/60480")
print("  D^8: (e2^4 - e2*delta^2)/1814400")
print("Cayley-block logarithmic derivative reduction checked exactly")
print("disc(x^mu-t) trace formula checked for 2 <= mu <= 9 and t in {2,3}")
print("persistent-fat-point and mixed-support discriminants vanish as claimed")
print("all exact checks passed")