"""Independent exact checks for T556, Yablonskii–Vorob’ev polynomials.

    sage -python checks.py

Requires generate.py and the campaign's check.py in the same directory.
Reads values.json when present (exported stored entries), otherwise checks
values produced by generate.py. No entries are written. Exact reconstruction
uses Bertola and Bothner, arXiv:1401.1408, equations (1.3)-(1.4), a finite
exponential coefficient sum and a fraction-free determinant. The nonlinear
recurrence in generate.py is not used by the determinant reconstruction.
"""
import json
from pathlib import Path
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
from check import exactness, measure, agrees_with, names_its_rings

S = PolynomialRing(QQ, 'x')
x = S.gen()


def determinant(rows):
    """Bareiss elimination over QQ[x], with exact remainders and row swaps."""
    a = [[S(v) for v in row] for row in rows]
    n = len(a)
    if n == 0:
        return S.one()
    previous = S.one()
    sign = 1
    for k in range(n-1):
        pivot_row = next((i for i in range(k,n) if a[i][k]), None)
        if pivot_row is None:
            return S.zero()
        if pivot_row != k:
            a[k],a[pivot_row] = a[pivot_row],a[k]
            sign = -sign
        pivot = a[k][k]
        for i in range(k+1,n):
            for j in range(k+1,n):
                numerator = pivot*a[i][j]-a[i][k]*a[k][j]
                quotient,remainder = numerator.quo_rem(previous)
                assert not remainder, 'Inexact Bareiss division'
                a[i][j] = quotient
            a[i][k] = S.zero()
        previous = pivot
    return sign*a[-1][-1]


def coefficient(k):
    """Coefficient of t^k in exp(x*t-4*t^3/3); zero for negative k."""
    return sum((QQ(-4)**j / (QQ(3)**j*QQ(factorial(j))*QQ(factorial(k-3*j)))
                *x**(k-3*j) for j in range(k//3+1)), S.zero())


def independent(n):
    rows = [[coefficient(n-2*i+j) for j in range(n)] for i in range(n)]
    factor = 1
    for k in range(1,n+1):
        factor *= (2*k+1)**(n-k)
    return factor*determinant(rows)


def painleve_residual(previous,current,n):
    """Numerator of w''-2*w^3-x*w-n for w=(log(previous/current))'."""
    a = previous.derivative()*current-previous*current.derivative()
    b = previous*current
    c = a.derivative()*b-a*b.derivative()
    d = c.derivative()*b-2*c*b.derivative()
    return d-2*a**3-x*a*b**2-n*b**3


def controls():
    assert determinant([]) == 1
    assert determinant([[0,1],[1,0]]) == -1
    assert determinant([[1,2],[2,4]]) == 0
    assert determinant([[x,1,0],[0,x,1],[1,0,x]]) == x**3+1
    assert coefficient(-1) == 0
    assert coefficient(0) == 1
    assert coefficient(3) == x**3/QQ(6)-QQ(4)/QQ(3)
    assert independent(2) == x**3+4
    assert painleve_residual(S.one(),x,1) == 0  # w=-1/x
    assert painleve_residual(x,x**3+5,2) != 0  # deliberate wrong constant


def check_values(values):
    controls()
    assert sorted(values)==list(range(13)), 'Expected the complete segment 0..12'
    assert not exactness(values)
    assert not agrees_with(values, independent), 'Determinant disagrees'
    # Roffelsen (2010), Table 1; initial Q_0 and Q_1 from its definition.
    reference={0:S.one(),1:x,2:x**3+4,3:x**6+20*x**3-80,
          4:x**10+60*x**7+11200*x,
          5:x**15+140*x**12+2800*x**9+78400*x**6-3136000*x**3-6272000,
          6:x**21+280*x**18+18480*x**15+627200*x**12-17248000*x**9
            +1448832000*x**6+19317760000*x**3-38635520000,
          7:x*(x**27+504*x**24+75600*x**21+5174400*x**18+62092800*x**15
               +13039488000*x**12-828731904000*x**9-49723914240000*x**6-3093932441600000),
          8:x**36+840*x**33+240240*x**30+32771200*x**27+2018016000*x**24
            +124309785600*x**21-6629855232000*x**18+407736096768000*x**15
            +126696533483520000*x**12+1769729356595200000*x**9
            +37164316488499200000*x**6-743286329769984000000*x**3-991048439693312000000}
    for n,p in values.items():
        assert p.degree()==n*(n+1)//2
        assert p.leading_coefficient()==1
        assert all(c.denominator()==1 for c in p.coefficients())
        if n in reference:
            assert p==reference[n], 'Published reference differs at n=%s: computed=%s; transcribed=%s' % (n,p,reference[n])
        if n:
            assert painleve_residual(values[n-1],p,n)==0
        if 0<n<12:
            assert values[n+1]*values[n-1]==x*p*p-4*(p*p.derivative(2)-p.derivative()**2)
    wrong_q5 = values[5]+2822400*x**3
    q4=values[4]
    assert wrong_q5*values[3]-(x*q4*q4-4*(q4*q4.derivative(2)-q4.derivative()**2)) == 2822400*x**3*values[3]
    assert painleve_residual(values[4],wrong_q5,5) != 0
    assert values[0]==1
    return {'table':'T556','entries':len(values),'independent_determinants':13,
            'painleve_II_identities':12,'stored_triple_recurrences':11,
            'Roffelsen_initial_polynomials':9,'controls':'all passed, including wrong-input rejection',
            'exactness':'all coefficients exact integers','measurement':measure(values),
            'skill_manifest':'../20261027T182939Z-build-skills/manifest.json',
            'skill_sha256':'e37efb7e2ededf6c1d7134c673f152f648f26e0267c556aebdf3ec65765d8b4c'}


if __name__=='__main__':
    source=Path(__file__).with_name('values.json')
    if source.exists():
        values={int(n):S(p) for n,p in json.loads(source.read_text()).items()}
        mode='stored readback'
    else:
        from generate import YablonskiiVorobev
        g=YablonskiiVorobev()
        values={int(p['n']):S(g.value(p,100)) for p in g.enumerate()}
        mode='local prefill'
        assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
    report=check_values(values)
    report['mode']=mode
    print(json.dumps(report,indent=2))
    print('VALUES_JSON='+json.dumps({n:str(p) for n,p in values.items()}))