"""Independent checks for Steklov eigenvalues -- numberdb.org/T489.

Run: sage -python checks.py
For a private draft set NUMBERDB_API_KEY, then:
     sage -python checks.py --stored
The generator must be beside this file. In stored mode every decimal is read
from the API, parsed as its exact last-place interval and checked against an
independent 180-digit computation. No write or review action is performed.

Controls: disk u=x has Rayleigh quotient 1; on the unit square,
u=(x-1/2)(y-1/2) has quotient 2. Positive rectangle modes are recomputed using
mpmath's trigonometric Robin quotients at 140 and 180 digits, independently
from the generator's pole-free Arb determinants. A complete separate branch
enumeration is sorted before comparison. Energy integrals use analytic
antiderivatives of squared factors and derivatives. Annulus eigenvalues are
independently obtained from the weighted boundary Dirichlet-to-Neumann matrix,
then tested on r^m + B r^-m with the inner normal -d/dr.
The square separation parameters are also checked against the beam roots in
T486. Its saved 100-digit entries are compared at their actual last-place
intervals; no assumption about equality of rounded midpoints is made.
"""
import json, os, sys
from decimal import Decimal
import mpmath as mp
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
import generate as gen


def decimal_rational(s):
    d=Decimal(str(s));t=d.as_tuple()
    n=int(''.join(map(str,t.digits)) or '0')*(-1 if t.sign else 1)
    return QQ(n)*QQ(10)**t.exponent


def written_interval(s):
    s=str(s)
    if '.' not in s and 'e' not in s.lower():
        q=QQ(s);return q,q
    q=decimal_rational(s);ulp=QQ(10)**Decimal(s).as_tuple().exponent
    return q-ulp,q+ulp


def inside_written(s,x):
    lo,hi=written_interval(s)
    q=decimal_rational(mp.nstr(x,175))
    assert lo<=q<=hi,(s,mp.nstr(x,110))


def mp_roots(a,b,dps):
    with mp.workdps(dps):
        out=[(mp.mpf(0),None)]
        if a==b:out.append((1/a,None))
        for ta,hb in [(a,b),(b,a)]:
            for trig in ('sin','cos'):
                for hyp in ('cosh','sinh'):
                    for k in range(6):
                        if trig=='sin' and hyp=='sinh' and k==0 and hb<=ta:continue
                        left=k*mp.pi/ta if trig=='sin' else (k+mp.mpf('.5'))*mp.pi/ta
                        right=left+mp.pi/(2*ta)
                        def f(x):
                            t=x*mp.cot(ta*x) if trig=='sin' else -x*mp.tan(ta*x)
                            h=x*mp.tanh(hb*x) if hyp=='cosh' else x*mp.coth(hb*x)
                            return t-h
                        # These quotients decrease, with opposite endpoint signs.
                        lo=left+mp.mpf('1e-6');hi=right-mp.mpf('1e-6')
                        assert f(lo)>0 and f(hi)<0
                        # mpmath Anderson/bisection is independent of the Arb code.
                        alpha=mp.findroot(f,(lo,hi),solver='anderson',tol=mp.mpf(10)**(-dps+15),maxsteps=1200)
                        assert left<alpha<right
                        sigma=alpha*(mp.tanh(hb*alpha) if hyp=='cosh' else mp.coth(hb*alpha))
                        out.append((+sigma,(+alpha,ta,hb,trig,hyp,k)))
        out.sort(key=lambda row:row[0])
        return out[:20]


def squares_integrals(alpha,length,kind):
    # Integral of F^2 and (F')^2 on [-length,length].
    x=alpha*length
    if kind=='cos':return length+mp.sin(2*x)/(2*alpha),alpha**2*(length-mp.sin(2*x)/(2*alpha))
    if kind=='sin':return length-mp.sin(2*x)/(2*alpha),alpha**2*(length+mp.sin(2*x)/(2*alpha))
    if kind=='cosh':return length+mp.sinh(2*x)/(2*alpha),alpha**2*(-length+mp.sinh(2*x)/(2*alpha))
    return -length+mp.sinh(2*x)/(2*alpha),alpha**2*(length+mp.sinh(2*x)/(2*alpha))


def quotient(alpha,a,b,trig,hyp):
    t,dt=squares_integrals(alpha,a,trig);h,dh=squares_integrals(alpha,b,hyp)
    tv=getattr(mp,trig)(alpha*a);hv=getattr(mp,hyp)(alpha*b)
    return (dt*h+t*dh)/(2*tv**2*h+2*hv**2*t)


def annulus_independent(dps):
    with mp.workdps(dps):
        result=[(mp.mpf(0),None),(mp.mpf(3)/(2*mp.log(2)),('radial',))]
        for m in range(1,25):
            # Boundary values at radii 1 and 2 give the Dirichlet-to-Neumann
            # matrix [[m*coth(mL), -m*csch(mL)],
            #         [-m*csch(mL)/2, m*coth(mL)/2]].
            # Multiplying boundary values by sqrt(radius) symmetrizes it.
            l=m*mp.log(2)
            matrix=mp.matrix([[m*mp.coth(l),-m*mp.csch(l)/mp.sqrt(2)],
                              [-m*mp.csch(l)/mp.sqrt(2),m*mp.coth(l)/2]])
            vals=mp.eigsy(matrix,eigvals_only=True)
            for j in range(2):result.extend([(+vals[j],(m,j))]*2)
        result.sort(key=lambda row:row[0]);return result[:20]


def main():
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
    stored='--stored' in sys.argv or os.environ.get('NUMBERDB_CHECK_STORED')=='1'
    # The wrapper forwards only selected flags; a tiny runner may call main.
    table=numberdb.table(gen.TABLE) if stored else None
    mp.mp.dps=180
    # Known nonzero controls are performed before the method is trusted.
    disk_control=mp.quad(lambda t:mp.cos(t)**2,[0,mp.pi,2*mp.pi])
    assert abs(mp.pi/disk_control-1)<mp.mpf('1e-170')
    a=mp.mpf('.5');area_energy=2*(2*a)*(2*a**3/3)
    boundary_mass=4*a*a*(2*a**3/3)
    assert area_energy/boundary_mass==2
    print('Controls: disk u=x gives 1; square bilinear mode gives 2.')
    independently={}
    for domain,a,b in [('unit-square',mp.mpf('.5'),mp.mpf('.5')),('rectangle-1x2',mp.mpf('.5'),mp.mpf(1))]:
        low=mp_roots(a,b,140);high=mp_roots(a,b,180)
        for (x,_),(y,mode) in zip(low,high):
            assert abs(x-y)<mp.mpf('1e-120')
            if mode:
                alpha,ta,hb,trig,hyp,k=mode
                assert abs(quotient(alpha,ta,hb,trig,hyp)-y)<mp.mpf('1e-155')
        independently[domain]=high
    for domain,a,b in [('unit-square',mp.mpf(1),mp.mpf(1)),('rectangle-1x2',mp.mpf(1),mp.mpf(2))]:
        scaled=mp_roots(a,b,160)
        for (x,_),(y,_) in zip(independently[domain],scaled):
            assert abs(x-2*y)<mp.mpf('1e-140')
    print('Scaling independently checked on all 40 rectangle modes.')
    alpha=mp.findroot(lambda x:mp.cos(x)*mp.cosh(x)+1,(1,2))
    assert abs(alpha*mp.tanh(alpha/2)-independently['unit-square'][1][0])<mp.mpf('1e-160')
    print('Displayed mpmath snippet checked.')
    independently['unit-disk']=[(mp.mpf(n//2),None) for n in range(1,21)]
    low=annulus_independent(140);high=annulus_independent(180)
    for (x,_),(y,mode) in zip(low,high):
        assert abs(x-y)<mp.mpf('1e-120')
        if mode and mode[0]!='radial':
            m,j=mode
            beta=(m+y)/(m-y)
            # Inner boundary: -u'(1)=sigma*u(1).
            assert abs((-m+m*beta)-y*(1+beta))<mp.mpf('1e-155')
            # Outer boundary: +u'(2)=sigma*u(2).
            residual=m*(mp.mpf(2)**(m-1)-beta*mp.mpf(2)**(-m-1))-y*(mp.mpf(2)**m+beta*mp.mpf(2)**(-m))
            assert abs(residual)<mp.mpf('1e-145')
    radial=mp.mpf(3)/(2*mp.log(2))
    # u=log(r)-1/sigma meets both radial boundary conditions.
    assert abs(mp.mpf('.5')-radial*(mp.log(2)-1/radial))<mp.mpf('1e-170')
    independently['annulus-1-2']=high
    minimum_digits=10000;count=0;written={}
    for domain in gen.DOMAINS:
        rows=gen.spectrum(domain)
        for n,(row,(x,_)) in enumerate(zip(rows,independently[domain]),1):
            if stored:
                r=table['Numbers'][domain][str(n)]
                text=r['number'] if isinstance(r,dict) else r
            else:text=to_text(row['number'],100)
            inside_written(text,x);count+=1;written[domain+','+str(n)]=str(text)
            if hasattr(row['number'],'rad') and row['number'].rad():
                v=row['number'];assert v.is_finite()
                rel=QQ(v.rad().exact_rational())/abs(QQ(v.mid().exact_rational()))
                digits=-mp.log10(mp.mpf(str(rel.numerator()))/int(rel.denominator()))
                minimum_digits=min(minimum_digits,float(digits))
                assert digits>110
        print(domain,'first values:',[mp.nstr(v,15) for v,_ in independently[domain][:6]])
    # External family comparison: read the actual beam table, not its generator.
    beams=numberdb.table('T486')['Numbers']
    compared=0
    for row in gen.spectrum('unit-square'):
        if 'alpha' not in row:continue
        alpha=row['alpha']; found=0
        for boundary in ['clamped-free','clamped-clamped']:
            for r in beams[boundary].values():
                s=r['number'] if isinstance(r,dict) else r
                lo,hi=written_interval(s)
                B=alpha.parent()
                if alpha.overlaps(B(lo).union(B(hi))):found+=1
        assert found==1
        compared+=1
    print('Independent decimal interval checks:',count)
    print('Square separation parameters matched to T486:',compared)
    print('Minimum relative enclosure digits:',round(minimum_digits,3))
    print('Stored mode:',stored)
    print('All controls, precision comparisons, energy and boundary checks passed.')
    print('WRITTEN_JSON',json.dumps(written,sort_keys=True))

if __name__=='__main__':main()