"""T712 independent check: periodic sequence sums and interval LDL^T.

Run with SageMath and numberdb: sage -python independent_check.py INPUT.json
INPUT is offline-results.json (builder evidence) or stored-input.json made
from the live table's stored entries. Without an argument use stored-input.json
if present, otherwise offline-results.json. Output consists of durable JSON
checkpoints and RESULT_JSON; save stdout, and extract RESULT_JSON to
check-results.json. No credentials, network, or writes to the website.

For period L, W_L(r)=(psi((r+1)/L)-psi(r/L))/L sums the exact weights
1/[m(m+1)] over m=r mod L. Direct periodic fractional products produce the
Gram and target. This does not use the generator's autocorrelation formula,
cotangent evaluations, or its linear solver. Exact rational arguments are
evaluated by complex Arb psi at 640 bits; values are real on these arguments.
Positive interval LDL pivots certify each finite solve. The two methods share
Sage/Arb arithmetic and the mathematical Hilbert-space definition.

Resume: extract CHECKPOINT_JSON lines into a list periodic-checkpoint.json,
place beside this script; completed Gram entries are read as outward ball
strings. These checkpoints are not themselves independent verdicts.
See skill-manifest.json, matching bundle
 d6661c2f21613524889fa3a57d3c95f8330344ca62e5a856211c08a282b7dbab.
"""
import os,sys,json,math,time,functools,hashlib
from pathlib import Path
from fractions import Fraction
from decimal import Decimal
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
R=RealBallField(640); C=ComplexBallField(640)
BASE=Path(__file__).resolve().parent


def narrow(x):
    return x.is_finite() and R(x.rad()) < R(QQ(1)/10**110)


def contains_written(x, s):
    """Require the whole independent ball inside the exact stored interval."""
    if not narrow(x): return False
    d=Decimal(s);mid=QQ(Fraction(s));unit=QQ(10)**d.as_tuple().exponent
    return bool(x-R(mid-unit)>0 and R(mid+unit)-x>0)


@functools.lru_cache(maxsize=None)
def weights(L):
    psi=[C(QQ(r)/L).psi().real() for r in range(1,L+2)]
    W=[(psi[r]-psi[r-1])/L for r in range(1,L+1)]
    assert all(narrow(x) and x>0 for x in W)
    assert (sum(W,R(0))-1).contains_zero()
    return W


def periodic_inner(i,j):
    L=math.lcm(i,j)
    return sum((R(QQ((r % i)*(r % j))/(i*j))*wr
                for r,wr in enumerate(weights(L),1)),R(0))


def target(k):
    return sum((R(QQ(r % k)/k)*wr for r,wr in enumerate(weights(k),1)),R(0))


def squarefree(k):
    # Independent elementary enumeration, rather than the generator's Sage test.
    return all(k % (d*d) for d in range(2,math.isqrt(k)+1))


def solve_ldl(A,b):
    n=len(b);L=[[R(i==j) for j in range(n)] for i in range(n)];D=[]
    for j in range(n):
        pivot=A[j][j]-sum((L[j][k]**2*D[k] for k in range(j)),R(0))
        assert narrow(pivot) and pivot>0
        D.append(pivot)
        for i in range(j+1,n):
            L[i][j]=(A[i][j]-sum((L[i][k]*L[j][k]*D[k] for k in range(j)),R(0)))/D[j]
    y=[]
    for i in range(n):y.append(b[i]-sum((L[i][j]*y[j] for j in range(i)),R(0)))
    z=[y[i]/D[i] for i in range(n)];x=[R(0)]*n
    for i in range(n-1,-1,-1):x[i]=z[i]-sum((L[j][i]*x[j] for j in range(i+1,n)),R(0))
    return x,D


def main():
    start=time.monotonic()
    filename=Path(sys.argv[1]) if len(sys.argv)>1 else BASE/('stored-input.json' if (BASE/'stored-input.json').exists() else 'offline-results.json')
    data=json.loads(filename.read_text());rows=data['entries'];nmax=max(int(r['params']['N']) for r in rows)
    expected=[{'support':s,'N':n} for s in ['all','squarefree'] for n in range(2,nmax+1)]
    assert [r['params'] for r in rows]==expected
    assert [k for k in range(2,31) if squarefree(k)]==[2,3,5,6,7,10,11,13,14,15,17,19,21,22,23,26,29,30]
    assert narrow(R(1)) and (sum(weights(1),R(0))-1).contains_zero()
    K22=periodic_inner(2,2);b2=target(2);control=1-b2*b2/K22
    assert K22.overlaps(R(2).log()/4) and b2.overlaps(R(2).log()/2)
    assert narrow(control) and control.overlaps(1-R(2).log())
    small=next(r['number'] for r in rows if r['params']=={'support':'all','N':2})
    assert contains_written(control,small)
    assert not contains_written(control+R(QQ(1)/10**6),small)
    assert not contains_written(R('nan'),small)
    assert not contains_written(R('0 +/- 2'),small)
    checkpoint=BASE/'periodic-checkpoint.json';memo={}
    if checkpoint.exists():
        for obj in json.loads(checkpoint.read_text()):
            assert obj['precision_bits']==640 and obj['method']=='periodic-digamma-v1'
            for k,v in obj['pairs'].items():memo[tuple(map(int,k.split(',')))]=R(v)
    for i in range(2,nmax+1):
        for j in range(i,nmax+1):
            if (i,j) not in memo:memo[i,j]=periodic_inner(i,j)
        print('CHECKPOINT_JSON '+json.dumps({'method':'periodic-digamma-v1','precision_bits':640,'i':i,'pairs':{f'{i},{j}':str(memo[i,j]) for j in range(i,nmax+1)}}),flush=True)
    b={k:target(k) for k in range(2,nmax+1)}
    assert all(narrow(x) and x.overlaps(R(k).log()/k) for k,x in b.items())
    gram_checks=0
    if 'gram' in data:
        for i in range(2,nmax+1):
            for j in range(i,nmax+1):
                actual=memo[i,j];generated=R(data['gram'][i-2][j-2])
                assert narrow(actual) and actual.overlaps(generated)
                gram_checks+=1
    done={}; comparisons=[];failures=[];minp=R(1);maxrad=R(0);minacc=99999
    for row in rows:
        p=row['params'];N=int(p['N']);s=p['support'];S=tuple(k for k in range(2,N+1) if s=='all' or squarefree(k))
        if S not in done:
            A=[[memo[min(i,j),max(i,j)] for j in S] for i in S];v=[b[k] for k in S]
            a,pivots=solve_ldl(A,v)
            x=1-sum((v[i]*a[i] for i in range(len(S))),R(0))
            assert narrow(x) and x>0 and x<1
            assert all((sum((A[i][j]*a[j] for j in range(len(S))),R(0))-v[i]).contains_zero() for i in range(len(S)))
            if 'systems' in data:
                old=data['systems'][','.join(map(str,S))]['optimizer']
                assert all(narrow(z) and z.overlaps(R(t)) for z,t in zip(a,old))
            done[S]=(x,a,pivots)
        x,a,pivots=done[S]
        ok=contains_written(x,row['number'])
        if not ok:failures.append({'params':p,'reason':'independent enclosure outside stored interval'})
        comparisons.append({'params':p,'number':row['number'],'check_ball':str(x),'contained':ok,'accuracy_bits':int(x.accuracy())})
        minp=min(minp,*pivots);maxrad=max(maxrad,R(x.rad()));minacc=min(minacc,int(x.accuracy()))
    # Nested-space inequalities, exact support plateaus, and inter-support bound.
    by={(r['params']['support'],r['params']['N']):R(r['check_ball']) for r in comparisons}
    for s in ['all','squarefree']:
        for N in range(3,nmax+1):
            if s=='squarefree' and not squarefree(N):assert by[s,N].overlaps(by[s,N-1])
            else:assert by[s,N-1]-by[s,N]>0
    for N in range(2,nmax+1):
        if N<=3:assert by['all',N].overlaps(by['squarefree',N])
        else:assert by['squarefree',N]-by['all',N]>0
    result={'table':'T712','input_sha256':hashlib.sha256(filename.read_bytes()).hexdigest(),
      'skill_manifest_sha256':hashlib.sha256((BASE/'skill-manifest.json').read_bytes()).hexdigest(),
      'checked_entries':len(rows),'total_entries':len(rows),'failures':failures,
      'known_answer_control':True,'negative_control_rejected':True,
      'controls':{'empty_span':1,'N2':str(control),'perturbation':'1e-6','nan_rejected':True,'vacuous_ball_rejected':True},
      'method':'periodic digamma weight sums; independent interval LDL transpose solve',
      'precision_bits':640,'minimum_pivot':str(minp),'maximum_radius':str(maxrad),
      'minimum_accuracy_bits':minacc,'gram_pairs_checked':gram_checks,
      'target_checks':nmax-1,'optimizer_systems_checked':len(done) if 'systems' in data else 0,
      'comparisons':comparisons,'elapsed_seconds':time.monotonic()-start}
    print('RESULT_JSON '+json.dumps(result),flush=True)
    if failures:raise SystemExit(1)

if __name__=='__main__':main()