verify_flavor_normalization.py

#!/usr/bin/env python3
"""Exact dependency-free checks for the flavor-normal Bethe root block.

Polynomials in d1,d2 are represented by coefficient dictionaries. The third
flavor coordinate is d3=-d1-d2. No numerical theta-function evaluation or
external computer-algebra package is used.
"""

from fractions import Fraction
from math import factorial


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


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


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(k: int) -> Poly:
    answer: Poly = {}
    for variable in deltas:
        answer = add(answer, power(variable, k))
    return answer


expected = {
    1: {},
    3: scale(delta, 3),
    5: scale(multiply(e2, delta), -5),
    7: scale(multiply(power(e2, 2), delta), 7),
    9: multiply(
        add(scale(power(e2, 3), -9), scale(power(delta, 2), 3)),
        delta,
    ),
}

for degree, target in expected.items():
    assert power_sum(degree) == target, degree

# The coefficient of ell^(k)(z) in H/delta is
# -2 * s_k / (k! * delta). The identities above imply the following exact
# coefficients, recorded as polynomials in e2 and delta.
ell_coefficients = {
    3: "-1",
    5: "e2/12",
    7: "-e2^2/360",
    9: "e2^3/20160 - delta^2/60480",
}
wp_coefficients = {
    1: "1",
    3: "-e2/12",
    5: "e2^2/360",
    7: "-e2^3/20160 + delta^2/60480",
}

# Independent rational checks of the displayed factorial reductions.
assert Fraction(-2 * 3, factorial(3)) == -1
assert Fraction(-2 * -5, factorial(5)) == Fraction(1, 12)
assert Fraction(-2 * 7, factorial(7)) == Fraction(-1, 360)
assert Fraction(-2 * -9, factorial(9)) == Fraction(1, 20160)
assert Fraction(-2 * 3, factorial(9)) == Fraction(-1, 60480)

print("constraint: d1 + d2 + d3 = 0")
print("e2 = d1*d2 + d2*d3 + d3*d1")
print("delta = d1*d2*d3")
print("s_1 = 0")
print("s_3 = 3*delta")
print("s_5 = -5*e2*delta")
print("s_7 = 7*e2^2*delta")
print("s_9 = (-9*e2^3 + 3*delta^2)*delta")
print("coefficients of ell^(k) in H/delta:")
for degree, coefficient in ell_coefficients.items():
    print(f"  k={degree}: {coefficient}")
print("normalized root block H/delta in Weierstrass derivatives:")
for derivative_order, coefficient in wp_coefficients.items():
    print(f"  wp^({derivative_order}): {coefficient}")
print("all exact identities passed")