"""Independent exact checks of Lamé spectral polynomials -- numberdb.org/T543.

Run: sage -python checks.py

Maier, arXiv:math-ph/0309005v4, Table 3, equations (8),(9),(53),
and Definition 7.1 supply an external coefficient-by-coefficient comparison.
The constants here are transcribed independently of generate.py. They use
Weierstrass B, not the generator's Jacobi h. Exact polynomial equality checks
all coefficients, not sampled moduli. Free-particle degenerations are checked
against h*product((h-j^2)^2), whose frequencies are the integer Fourier modes.
No approximate arithmetic or fitted data enter these checks.
"""
import json
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing

CHECK_RING = PolynomialRing(QQ, names=('h', 'm'), order='lex')
x, q = CHECK_RING.gens()


def maier_polynomial(n):
    """Translate independently tabulated Weierstrass factors to Jacobi form."""
    B = -x + QQ(n*(n+1))/3*(q+1)
    g2 = QQ(4)/3*(q*q-q+1)
    g3 = QQ(4)/27*(q-2)*(2*q-1)*(q+1)
    type_i = {
        1: CHECK_RING.one(),
        2: B**2-3*g2,
        3: B,
        4: B**3-52*g2*B+560*g3,
        5: B**2-27*g2,
        6: B**4-294*g2*B**2+7776*g3*B+3465*g2**2,
        7: B**3-196*g2*B+2288*g3,
        8: B**5-1044*g2*B**3+48816*g3*B**2+112320*g2**2*B-4665600*g2*g3,
    }
    result = -type_i[n]  # Degree 2n+1 is odd and B reverses the sign of h.
    for e in [(2-q)/3, (2*q-1)/3, -(q+1)/3]:
        type_ii = {
            1: B-e,
            2: B+3*e,
            3: B**2-6*e*B+45*e**2-15*g2,
            4: B**2+10*e*B-35*e**2-7*g2,
            5: B**3-15*e*B**2+(315*e**2-132*g2)*B+QQ(675)/4*e*g2+QQ(2835)/4*g3,
            6: B**3+21*e*B**2+(-189*e**2-84*g2)*B-QQ(3465)/4*e*g2+QQ(4455)/4*g3,
            7: B**4-28*e*B**3+(1134*e**2-574*g2)*B**2+(3409*e*g2+8525*g3)*B-QQ(292383)/4*e**2*g2-QQ(175175)/4*e*g3+22113*g2**2,
            8: B**4+36*e*B**3+(-594*e**2-414*g2)*B**2+(-9855*e*g2+12285*g3)*B+QQ(245025)/4*e**2*g2+QQ(552825)/4*e*g3+7425*g2**2,
        }
        result *= type_ii[n]
    return result


def check_values(values):
    # Control first: the external conversion must reproduce a known case.
    assert maier_polynomial(1) == (x-1)*(x-q)*(x-q-1)
    assert maier_polynomial(2) == (x*x-4*(q+1)*x+12*q)*(x-q-1)*(x-4*q-1)*(x-q-4)
    rows = []
    for key, value in values.items():
        n = int(key)
        p = CHECK_RING(value)
        N = n*(n+1)
        assert p == maier_polynomial(n), ('Maier coefficients', n)
        assert p.degree(x) == 2*n+1
        assert p.coefficient({x:2*n+1}) == 1
        assert all(v.denominator() == 1 for v in p.coefficients())
        free = x
        for j in range(1,n+1):
            free *= (x-j*j)**2
        assert p(x,0) == free, ('free-particle limit', n)
        assert p == -p(N-x,1-q), ('complementary symmetry',n)
        assert p(QQ(N)/2,QQ(1)/2) == 0, ('lemniscatic central root',n)
        # The m=1 degeneration follows from the independently known m=0 limit.
        degenerate = x-N
        for j in range(1,n+1):
            degenerate *= (x-N+j*j)**2
        assert p(x,1) == degenerate
        assert p.coefficient({x:2*n}) == -QQ(N*(2*n+1))/3*(1+q)
        rows.append({'n':n,'degree_h':int(p.degree(x)),
                     'terms':len(p.monomials()),'characters':len(str(p)),
                     'maier_equal':True,'identities_pass':True})
    return {'control_n1_n2':True,'all_pass':True,'rows':rows}


if __name__ == '__main__':
    from generate import LameSpectral
    generator = LameSpectral()
    values = {p['n']:generator.value(p) for p in generator.enumerate()}
    print(json.dumps(check_values(values),indent=2))