"""Independent rational-breakpoint check for the Baez-Duarte R_N errors.

sage -python independent_check.py stored.json [result.json]
Input is a NumberDB document or an offline mapping N -> encoded interval.
This does not import generate.py. It builds Mobius values using the divisor
identity sum_{d|n} mu(d) = [n=1], builds the periodic step function from exact
jumps, and integrates its square and logarithmic cross term by endpoint
values of digamma and the deflated Hurwitz-zeta expansion at s=1.
See derivation.md. Shared with the generator: Arb, exact rationals, and the
mathematical definition. Different arithmetic construction, kernel and cross
term implementation. This is builder evidence, not an independent actor verdict.
"""
import json,sys
from functools import lru_cache
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.integer_ring import ZZ
from sage.rings.real_arb import RealBallField
from sage.rings.real_mpfi import RealIntervalField
from sage.rings.complex_arb import ComplexBallField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from decimal import Decimal
BITS=500
R=RealBallField(BITS);C=ComplexBallField(BITS)
P=PolynomialRing(C,'s');s=P.gen()

@lru_cache(None)
def endpoint(a):
    z=(1+s)._zeta_series(2,a=C(a),deflate=True)
    assert z[0].imag().contains_zero() and z[1].imag().contains_zero()
    # z[0] = -psi(a); z[1] = -gamma_1(a)
    return -z[0].real(),z[1].real()

def weights(n):
    mu=[0]*(n+1)
    for j in range(1,n+1):
        mu[j]=(1 if j==1 else 0)-sum(mu[d] for d in range(1,j) if j%d==0)
    M=[sum(mu[:j+1]) for j in range(n+1)]
    w={k:QQ(M[n//k])/k for k in range(1,n) if M[n//k]}
    assert sum((k*v for k,v in w.items()),QQ(0))==0
    return w

def direct(n):
    n=int(n);w=weights(n);events={QQ(0):QQ(0),QQ(1):QQ(0)}
    for k,v in w.items():
        for j in range(1,k):
            a=QQ(j)/k;events[a]=events.get(a,QQ(0))-v
    points=sorted(events);c=QQ(0);norm=R(0);minus_cross=R(0)
    for a,b in zip(points,points[1:]):
        c+=events[a]
        if not c:continue
        assert a>=QQ(1)/n
        pa,za=endpoint(a);pb,zb=endpoint(b);delta=pb-pa
        norm+=R(c*c)*delta/n
        minus_cross+=R(c)*((R(n).log()+1)*delta+zb-za)/n
    return 2+norm+2*minus_cross

def control3():
    # F_3=-1/2 on [1/2,1); the digamma duplication identity gives
    # norm=log(2)/6 and the elementary cancelled Mellin cross term below.
    l2=R(2).log()
    return 2+l2/6-2*l2*(1-R.euler_constant()+R(3).log()-l2/2)/3

def compare(encoded,value):
    if not value.is_finite() or value.accuracy()<380:return False
    dec=Decimal(str(encoded).rstrip('?'))
    num,den=dec.as_integer_ratio(); center=QQ(num)/den
    ulp=QQ(10)**int(dec.as_tuple().exponent)
    return bool(value.lower()>=center-ulp and value.upper()<=center+ulp)

def run(input_path,result_path):
    doc=json.loads(Path(input_path).read_text())
    if 'Numbers' in doc:
        rows={str(k):v['number'] for k,v in doc['Numbers'].items()}
    else:rows=doc
    result={'checked_entries':0,'total_entries':len(rows),'failures':[],
            'known_answer_control':False,'negative_control_rejected':False,
            'method':'rational-breakpoint integration with deflated Hurwitz zeta, 500-bit Arb',
            'comparisons':[]}
    v3=direct(3)
    result['known_answer_control']=bool(direct(1)==2 and direct(2)==2 and (v3-control3()).contains_zero() and v3.accuracy()>380)
    result['negative_control_rejected']=not compare('1.573641014547909618621335725814',v3)
    result['invalid_ball_controls']=not compare('2.00',R('nan')) and not compare('2.00',R(RealIntervalField(BITS)(-100,100)))
    for n,encoded in rows.items():
        value=direct(int(n));ok=compare(encoded,value)
        result['checked_entries']+=1
        result['comparisons'].append({'N':int(n),'ok':ok,'ball':str(value),'accuracy_bits':int(value.accuracy())})
        if not ok:result['failures'].append(int(n))
        Path(result_path).write_text(json.dumps(result,indent=2)+'\n')
        if int(n)%10==0:print('checked through',n,flush=True)
    assert result['known_answer_control'] and result['negative_control_rejected'] and result['invalid_ball_controls'] and not result['failures'],result['failures']
    print('Independent check:',result['checked_entries'],'passed; controls passed')

if __name__=='__main__':run(sys.argv[1],sys.argv[2] if len(sys.argv)>2 else 'check-results.json')