"""Lamé band-edge eigenvalues -- numberdb.org/T544.

Run with SageMath (which includes NumPy and SciPy):
    sage -pip install numberdb              # once
    sage -python generate.py                 # verify every stored value
    sage -python generate.py --fill          # fill existing draft, API key in env

This file does not create, offer, review, or publish a table. The client method
named publish writes entries and attached files to the existing table only.

The Jacobi operator is -d²/dx²+n(n+1)m sn²(x|m), with m the parameter.
Exact rational gauge transformations construct all invariant polynomial sectors.
Their characteristic polynomials are computed by a tridiagonal determinant
recurrence. Rational factorization detects exact rational eigenvalues. All other
roots are isolated by exact Sturm counts and refined by rational bisection.
A nonzero-width MPFI interval is formed from the exact rational endpoints,
with outward rounding. No approximate root finder determines the answer.

The default selection and its finite cutoffs are explained in the attached
argument-selection.md. ROOT_GUARD_BITS=80 gives more than 120 decimal places
of absolute enclosure at digits=100; the checks measure the actual widths.
"""
import os
import sys
from functools import lru_cache
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
from sage.rings.real_mpfi import RealIntervalField

ROOT_GUARD_BITS = 80
R = PolynomialRing(QQ, 't')
t = R.gen()
S = PolynomialRing(QQ, 'h')
h = S.gen()


def sector_polynomials(n, m):
    """Characteristic polynomial of H on each admissible sn^a cn^b dn^c p(t)."""
    n, m = int(n), QQ(m)
    F = 4*t*(1-t)*(1-m*t)
    G = 2*(1-2*(1+m)*t+3*m*t*t)
    out = []
    for a,b,c in product((0,1), repeat=3):
        if n < a+b+c or (n-a-b-c) % 2:
            continue
        d = (n-a-b-c)//2
        q = QQ(a)/(2*t) - QQ(b)/(2*(1-t)) - QQ(c)*m/(2*(1-m*t))
        B = R(-2*F*q-G)
        C = R(n*(n+1)*m*t-F*(q.derivative()+q*q)-G*q)
        columns = []
        for k in range(d+1):
            p = t**k
            image = -F*p.derivative(2)+B*p.derivative()+C*p
            assert image.degree() <= d
            assert all(image[r] == 0 for r in range(d+1) if abs(r-k)>1)
            columns.append(image)
        prevprev, prev = S(1), h-columns[0][0]
        for k in range(1,d+1):
            cur = (h-columns[k][k])*prev-columns[k-1][k]*columns[k][k-1]*prevprev
            prevprev, prev = prev, cur
        out.append(((a,b,c),prev))
    assert sum(f.degree() for _,f in out) == 2*n+1
    return out


@lru_cache(maxsize=None)
def spectral_polynomial(n,m):
    p = S(1)
    for _,f in sector_polynomials(n,QQ(m)):
        p *= f
    assert p.is_monic() and p.gcd(p.derivative()).degree()==0
    return p


def sturm_sequence(f):
    seq=[f,f.derivative()]
    while seq[-1].degree()>0:
        rem = -(seq[-2] % seq[-1])
        assert rem
        # Positive rescaling preserves every Sturm variation.
        seq.append(rem/abs(rem.leading_coefficient()))
    return seq


def variations(seq,x):
    signs=[1 if v>0 else -1 for f in seq if (v:=f(x)) != 0]
    return sum(a!=b for a,b in zip(signs,signs[1:]))


def isolate_factor(f,bound,bits):
    """All real roots of an irreducible nonlinear rational polynomial."""
    seq=sturm_sequence(f)
    left,right=QQ(-1),QQ(bound)
    vl,vr=variations(seq,left),variations(seq,right)
    assert vl-vr==f.degree()
    todo=[(left,right,vl,vr)]
    brackets=[]
    while todo:
        lo,hi,vlo,vhi=todo.pop()
        count=vlo-vhi
        if count==0:
            continue
        if count==1:
            brackets.append((lo,hi))
            continue
        mid=(lo+hi)/2
        assert f(mid)!=0  # nonlinear irreducible factors have no rational root
        vmid=variations(seq,mid)
        todo.extend([(lo,mid,vlo,vmid),(mid,hi,vmid,vhi)])
    target=QQ(1)/QQ(2)**bits
    out=[]
    for lo,hi in brackets:
        flo=f(lo)
        assert flo*f(hi)<0
        while hi-lo>target:
            mid=(lo+hi)/2
            fm=f(mid)
            assert fm!=0
            if flo*fm<0:
                hi=mid
            else:
                lo,flo=mid,fm
        assert variations(seq,lo)-variations(seq,hi)==1
        out.append((lo,hi,f))
    return out


@lru_cache(maxsize=None)
def root_certificates(n,m,digits=100):
    n,m=int(n),QQ(m)
    assert n>=1 and 0<m<1
    bits=numberdb.bits(digits,losing=ROOT_GUARD_BITS)
    roots=[]
    p=spectral_polynomial(n,str(m))
    for f,multiplicity in p.factor():
        assert multiplicity==1
        if f.degree()==1:
            r=-f[0]/f[1]
            assert p(r)==0
            roots.append((r,r,f))
        else:
            roots.extend(isolate_factor(f,2*n*(n+1)+1,bits))
    roots.sort(key=lambda r:r[0])
    assert len(roots)==2*n+1
    assert all(a[1]<b[0] for a,b in zip(roots,roots[1:]))
    return tuple(roots)


def band_edges(n,m,digits=100):
    field=RealIntervalField(numberdb.bits(digits,losing=ROOT_GUARD_BITS+32))
    return [lo if lo==hi else field(lo,hi)
            for lo,hi,_ in root_certificates(int(n),str(QQ(m)),digits)]


class LameBandEdges(numberdb.Generator):
    table='T544'
    parameters=('n','m','j')
    type='R'
    digits=100
    rigour='proven'
    files=('generate.py','checks.py','check.py','argument-selection.md','numerical-evidence.json','skill-manifest.json')

    def enumerate(self,lemniscatic_order=12,reference_order=6):
        for n in range(1,max(lemniscatic_order,reference_order)+1):
            for m,limit in [('1/2',lemniscatic_order),('1/4',reference_order),('3/4',reference_order)]:
                if n<=limit:
                    for j in range(2*n+1):
                        yield {'n':n,'m':m,'j':j}

    def value(self,params,digits):
        n,j=int(params['n']),int(params['j'])
        assert 0<=j<=2*n
        return band_edges(n,params['m'],digits)[j]


if __name__=='__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        key=sys.stdin.read().strip()
        if '=' in key and key.split('=',1)[0].isidentifier():
            key=key.split('=',1)[1].strip().strip('\"\'')
        numberdb.configure(api_key=key)
    generator=LameBandEdges()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1':
        print(generator.publish(message='Certified finite-gap edges from exact sector polynomials'))
    else:
        report=generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)