"""Lamé spectral polynomials -- numberdb.org/T543.

Run with SageMath (the numberdb Python package is required):
    sage -pip install numberdb
    sage -python generate.py                 # verify all stored entries
    sage -python checks.py                   # independent exact checks

This file produces the monic Jacobi-form spectral polynomial in QQ[h,m].
The finite matrices act on sn(x|m)^a cn(x|m)^b dn(x|m)^c times polynomials
in t=sn(x|m)^2. Their tridiagonal determinants are evaluated by continuants,
using only exact additions and multiplications. No root approximation enters.
The coefficient of the kinetic term is 1; m is the squared modulus.
See argument-selection.md for the finite order bound and source locators.
An authorized writer can fill the private draft with
LameSpectral().publish(); this operation submits entries, not public review.
"""
import sys
from itertools import product
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing

R = PolynomialRing(QQ, names=('h', 'm'), order='lex')
h, m = R.gens()
MAX_ORDER = 4  # 843 characters; order 5 already needs 1462, order 8 needs 5039.


def spectral_polynomial(n):
    """Return P_n(h,m), monic in h, for a positive integer n."""
    n = int(n)
    if n < 1:
        raise ValueError('The Lamé order must be a positive integer.')
    result = R.one()
    dimension = 0
    for a, b, c in product((0, 1), repeat=3):
        s = a + b + c
        if n < s or (n - s) % 2:
            continue
        d = (n - s) // 2
        dimension += d + 1
        previous, current = R.zero(), R.one()
        for k in range(d + 1):
            diagonal = (4*k*(k*(1+m) + a*(1+m) + b + c*m)
                        + (a+b)**2 + m*(a+c)**2)
            lower = -2*k*(2*k+2*a-1)
            upper_previous = m*(n*(n+1)-(2*(k-1)+s)*(2*(k-1)+s+1))
            following = (h-diagonal)*current-lower*upper_previous*previous
            previous, current = current, following
        result *= current
    assert dimension == 2*n+1
    return result


class LameSpectral(numberdb.Generator):
    table = 'T543'
    parameters = ('n',)
    type = 'Q[]'
    rigour = 'exact'

    def enumerate(self, max_order=MAX_ORDER):
        for n in range(1, int(max_order)+1):
            yield {'n': str(n)}

    def value(self, params, digits=100):
        return spectral_polynomial(int(params['n']))


if __name__ == '__main__':
    report = LameSpectral().verify(sample=None)
    print(report)
    sys.exit(0 if report.ok else 1)