"""Independent checks for T631. Run: sage -python checks.py.

The scalar check uses SO(4) trigonometric character traces and separately
regularizes the factors m-1 and m+1. It shares no multiplicity or continuation
code with generate.py. Ray-Singer torsion then determines an independent
coexact log determinant. mpmath is corroboration, not the rigorous enclosure.
With NUMBERDB_CHECK_STORED=1, also read T631, run Generator.verify on all
entries, and apply independent checks to the actual decimal intervals.
"""
import os
import sys
import json
from pathlib import Path
from functools import lru_cache
from math import gcd
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
from sage.rings.complex_arb import ComplexBallField
import mpmath as mp
from generate import CoexactLensDeterminants, cases, presentations, coefficients, logdet
from check import exactness, measure, names_its_rings

mp.mp.dps=175
HERE=Path(__file__).resolve().parent


def mpq(a,b=1):
    return mp.mpf(int(a))/int(b)


@lru_cache(None)
def hz(j,n,p):
    a=mpq(n,p)
    return mp.zeta(-j,a),mp.zeta(-j,a,derivative=1)


@lru_cache(None)
def scalar_traces(p,q):
    # Coefficients of m^2,m,1 on m = r mod p. The central antipodal
    # element contributes (-1)^(m-1)*m^2. Noncentral equal-angle
    # elements contribute m*sin(m*alpha)/sin(alpha).
    out=[]
    for j in range(p):
        rows=[]
        v=(q*j)%p
        alpha=2*mp.pi*j/p; beta=2*mp.pi*v/p
        for r in range(2,p+2):
            if j==0:
                row=(mp.mpf(1),mp.mpf(0),mp.mpf(0))
            elif 2*j==p:
                row=(mp.mpf((-1)**(r-1)),mp.mpf(0),mp.mpf(0))
            elif v==j or v==p-j:
                row=(mp.mpf(0),mp.sin(r*alpha)/mp.sin(alpha),mp.mpf(0))
            else:
                row=(mp.mpf(0),mp.mpf(0),(mp.cos(r*beta)-mp.cos(r*alpha))/(mp.cos(beta)-mp.cos(alpha)))
            rows.append(row)
        out.append(rows)
    return out


@lru_cache(None)
def scalar_logdet(p,q,u):
    traces=scalar_traces(p,q)
    total=mp.mpf(0); lp=mp.log(p); mean_linear=mp.mpf(0)
    for ri,r in enumerate(range(2,p+2)):
        A,B,C=[mp.fsum(mp.cos(2*mp.pi*u*j/p)*traces[j][ri][k] for j in range(p))/p for k in range(3)]
        mean_linear+=B
        for c in (-1,1):
            shifted=(C-c*B+c*c*A,B-2*c*A,A)
            for k in range(3):
                z,dz=hz(k,r+c,p)
                total+=shifted[k]*p**k*(lp*z-dz)
    assert abs(mean_linear)<mp.mpf('1e-165')
    return total


def torsion_log(p,q,u):
    if u==0:
        return mp.log(2*mp.pi**2/p**2)
    v=u*pow(q,-1,p)%p
    return mp.log(4*mp.sin(mp.pi*u/p)*mp.sin(mp.pi*v/p))


def rational_text(text):
    if '/' in text:
        a,b=text.split('/'); return mpq(a,b)
    return mp.mpf(text)


def scalar_multiplicity_direct(p,q,u,m):
    # Harmonic homogeneous polynomials: monomials of degree m-1 in
    # z1,z1bar,z2,z2bar, minus |z|^2 times degree m-3.
    def monomials(n):
        if n<0:return 0
        return sum((a-b+q*(c-(n-a-b-c))-u)%p==0
                   for a in range(n+1) for b in range(n-a+1)
                   for c in range(n-a-b+1))
    return monomials(m-1)-monomials(m-3)


def main():
    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('\"\'')
        os.environ['NUMBERDB_API_KEY']=key
    g=CoexactLensDeterminants(); R=RealBallField(700)
    counts={'presentations':len(list(presentations())),'space_character_pairs':len(list(cases()))}
    assert counts=={'presentations':46,'space_character_pairs':332}
    assert not names_its_rings(str(HERE/'generate.py'))
    # Known nonzero controls, before any broad comparison.
    expected=-3*R(3).zeta()/R.pi()**2-2*R(2).log()
    assert logdet(2,1,1).overlaps(expected)
    scalar_sphere=mp.log(mp.pi)+mp.zeta(3)/(2*mp.pi**2)
    assert abs(scalar_logdet(2,1,1)-(mp.log(2)-3*mp.zeta(3)/(2*mp.pi**2)))<mp.mpf('1e-160')
    assert abs(scalar_logdet(2,1,0)+scalar_logdet(2,1,1)-scalar_sphere)<mp.mpf('1e-160')
    print('Known coexact and scalar controls pass.',flush=True)
    n_sym=0; n_cover=0
    # Equality of all polynomial coefficients proves equality of every
    # multiplicity, not just a finite range of modes.
    for p in range(2,21):
        for q in range(1,p):
            if gcd(p,q)!=1:continue
            r=pow(q,-1,p)
            for u in range(p):
                c=coefficients(p,q,u)
                assert c==coefficients(p,q,(-u)%p)
                assert c==coefficients(p,p-q,u)
                assert c==coefficients(p,r,(u*r)%p)
                n_sym+=1
            allc=[coefficients(p,q,u) for u in range(p)]
            for i in range(2*p):
                m=2+i
                assert sum(v[i][1] for v in allc)==2*(m*m-1)
                assert sum(v[i][2] for v in allc)==8*m*p
                assert sum(v[i][3] for v in allc)==8*p*p
                n_cover+=1
    print('Exact symmetry and covering multiplicity identities pass.',flush=True)
    # Independent low-mode scalar monomial controls on every selected case.
    n_modes=0
    for p,q,u in cases():
        tr=scalar_traces(p,q)
        for m in range(2,8):
            ri=(m-2)%p
            value=mp.fsum(mp.cos(2*mp.pi*u*j/p)*(tr[j][ri][0]*m*m+tr[j][ri][1]*m+tr[j][ri][2]) for j in range(p))/p
            assert abs(value-scalar_multiplicity_direct(p,q,u,m))<mp.mpf('1e-160')
            n_modes+=1
    print('Independent scalar monomial controls pass.',flush=True)
    values={}; max_difference=mp.mpf(0); min_acc=100000; min_high_acc=100000
    results={}; scalar_results={}
    for p,q,u in cases():
        L=logdet(p,q,u,100); H=logdet(p,q,u,150)
        assert L.is_finite() and H.is_finite() and L.overlaps(H)
        other=2*(scalar_logdet(p,q,u)-torsion_log(p,q,u))
        diff=abs(mp.mpf(to_text(H,160))-other)
        max_difference=max(max_difference,diff)
        assert diff<mp.mpf('1e-145'),(p,q,u,mp.nstr(diff))
        key=f'{p},{q},{u}'
        results[key]={'independent_logdet':mp.nstr(other,155),'difference':mp.nstr(diff,6)}
        scalar_results[key]=mp.nstr(scalar_logdet(p,q,u),155)
        for convention,x,y in [('det',L.exp(),H.exp()),('logdet',L,H)]:
            assert x.is_finite() and y.is_finite()
            min_acc=min(min_acc,int(x.accuracy()));min_high_acc=min(min_high_acc,int(y.accuracy()))
            values[key+','+convention]=x
            assert x.overlaps(y)
        if u==0:print('checked',p,q,flush=True)
    assert not exactness(values)
    # Covering determinants from stored sectors, not used to generate them.
    cover_log=R(3).zeta()/R.pi()**2-2*(2*R.pi()).log()
    for p,q in presentations():
        total=R(0)
        for u in range(p):total+=logdet(p,q,min(u,p-u),100)
        assert total.overlaps(cover_log)
    # Read the existing published companions; archived documents are
    # optional for offline reproducibility, and never used in generation.
    companions={}
    for tid in ['T599','T595','T600']:
        path=HERE/(tid+'.json')
        companions[tid]=json.loads(path.read_text()) if path.exists() else numberdb.table(tid)
    for p,q,u in cases():
        if u:
            text=companions['T599']['Numbers'][str(p)][str(q)][str(u)]
            if isinstance(text,dict):text=text['number']
            assert abs(rational_text(text)-mp.exp(torsion_log(p,q,u)))<mp.mpf('1e-95')
    for tid,x in [('T595',scalar_sphere),('T600',scalar_logdet(2,1,0))]:
        text=companions[tid]['Numbers']['3']['logdet']
        if isinstance(text,dict):text=text['number']
        assert abs(rational_text(text)-x)<mp.mpf('1e-95')
    evidence={'table':'T631','skill_manifest':'../20261030T204529Z-build-skills/manifest.json',
              'skill_sha256':'10c1d8bd62d6ddfa2c016375ad82879151f0ab7a9516c2f371432108a85be009',
              'counts':counts,'entries':len(values),'working_bits':numberdb.bits(100,losing=180),
              'minimum_accuracy_bits':min_acc,'higher_precision_minimum_accuracy_bits':min_high_acc,
              'independent_mpmath_dps':mp.mp.dps,'maximum_logdet_discrepancy':mp.nstr(max_difference,12),
              'symmetry_cases':n_sym,'covering_polynomial_checks':n_cover,'scalar_monomial_checks':n_modes,
              'size':measure({k:to_text(v,100) for k,v in values.items()}),'results':results,'independent_scalar_logdets':scalar_results,
              'versions':g.environment(),'stored_checked':False}
    if os.environ.get('NUMBERDB_CHECK_STORED')=='1':
        report=g.verify(sample=None);print(report);assert report.ok
        document=numberdb.table('T631')
        n=0
        for p,q,u in cases():
            row=document['Numbers'][str(p)][str(q)][str(u)]
            other=mp.mpf(results[f'{p},{q},{u}']['independent_logdet'])
            for convention,x in [('det',mp.exp(other)),('logdet',other)]:
                text=row[convention]
                if isinstance(text,dict):text=text['number']
                from decimal import Decimal
                dec=Decimal(text); ulp=mp.power(10,dec.as_tuple().exponent)
                assert abs(mp.mpf(text)-x)<ulp/2,(p,q,u,convention)
                # A 10-ulp perturbation must fail the same acceptance rule.
                assert abs(mp.mpf(text)+10*ulp-x)>ulp
                n+=1
        evidence['stored_checked']=True;evidence['stored_entries_checked']=n
    print('EVIDENCE_JSON_START',flush=True)
    print(json.dumps(evidence,indent=2),flush=True)
    print('EVIDENCE_JSON_END',flush=True)


if __name__=='__main__':main()