"""Generalized Hermite polynomials for Painlevé IV -- numberdb.org/T557.

Run with SageMath (or a suitable passagemath installation):
    sage -pip install numberdb
    sage -python generate.py                    # verify every stored entry
    sage -python generate.py --fill             # write entries to this draft
    sage -python checks.py                      # independent exact checks

Set NUMBERDB_API_KEY for authenticated reads and writes of the private draft.
For the campaign wrapper, NUMBERDB_KEY_FROM_STDIN=1 reads the key from stdin;
NUMBERDB_PUBLISH=1 writes entries, but never publishes or offers the table.

Every value is a monic polynomial over QQ in the unchanged variable x.
Compute physicists' Hermite polynomials by their three-term recurrence, then
compute the Wronskian by fraction-free elimination with exact remainder tests.
No floating-point computation, numerical differentiation, symmetry reduction,
or Painlevé differential recurrence is used in this generator.
See argument-selection.md for the two discrete indices and finite degree bound.
"""
import os
import sys
from functools import lru_cache
from math import factorial

import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing

R = PolynomialRing(QQ, 'x')
x = R.gen()
MAX_DEGREE = 50


@lru_cache(None)
def hermite(k):
    if k == 0:
        return R.one()
    if k == 1:
        return 2*x
    return 2*x*hermite(k-1) - 2*(k-1)*hermite(k-2)


def determinant(rows):
    """Bareiss determinant in QQ[x], checking every exact division."""
    a = [list(row) for row in rows]
    size = len(a)
    if size == 0:
        return R.one()
    previous = R.one()
    sign = 1
    for k in range(size-1):
        if a[k][k] == 0:
            pivot_row = next((i for i in range(k+1,size) if a[i][k]), None)
            if pivot_row is None:
                return R.zero()
            a[k],a[pivot_row] = a[pivot_row],a[k]
            sign = -sign
        pivot = a[k][k]
        for i in range(k+1,size):
            for j in range(k+1,size):
                numerator = pivot*a[i][j] - a[i][k]*a[k][j]
                a[i][j], remainder = numerator.quo_rem(previous)
                assert remainder == 0, ('inexact Bareiss division',k,i,j)
            a[i][k] = R.zero()
        previous = pivot
    return sign*a[-1][-1]


@lru_cache(None)
def generalized_hermite(m,n):
    if m < 0 or n < 0:
        raise ValueError('Indices must be nonnegative')
    rows = [[hermite(m+j).derivative(i) for j in range(n)] for i in range(n)]
    polynomial = determinant(rows)
    leading = QQ(2)**(m*n+n*(n-1)//2)
    for j in range(1,n):
        leading *= factorial(j)
    assert polynomial.degree() == m*n
    assert polynomial.leading_coefficient() == leading
    return polynomial * (QQ(1)/leading)


class GeneralizedHermite(numberdb.Generator):
    table = 'T557'
    parameters = ('m','n')
    type = 'Q[]'
    rigour = 'exact'
    digits = 100  # Exact polynomial serialization does not round coefficients.
    files = ('generate.py','checks.py','check.py','argument-selection.md',
             'source-notes.md','skill-evidence.json','numerical-evidence.json')

    def enumerate(self, max_degree=MAX_DEGREE):
        for degree in range(1,max_degree+1):
            for m in range(1,degree+1):
                if degree % m == 0:
                    yield {'m':m,'n':degree//m}

    def value(self, params, digits):
        return generalized_hermite(int(params['m']),int(params['n']))


def authenticate():
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        key = sys.stdin.read().strip()
        if '=' in key and key.split('=',1)[0].replace('_','').isalnum():
            key = key.split('=',1)[1].strip().strip('\"\'')
        numberdb.configure(api_key=key)


if __name__ == '__main__':
    authenticate()
    generator = GeneralizedHermite()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Fill monic generalized Hermite polynomials through degree 50, with exact independent checks.'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)