"""Independent numerical checks for Values of the Riesz function, T645.

Run: sage -python checks.py [--stored]
Keep generate.py, check.py, skill-manifest.json and T644.json beside this file.
The optional --stored uses the authenticated NumberDB client to read T645;
otherwise all results are checked before any write. Output is JSON evidence.
No website mutations are performed by this script.

The independent method splits the Mobius representation at N=16. It evaluates
n<=N as exponentials and expands the rest through k=99. The residual coefficient
is 1/zeta(2k+2)-sum_{n<=N}mu(n)/n^(2k+2). Taylor's integral remainder for e^-t
on t>=0 is bounded by t^K/K!, hence the omitted tail of R(x) is at most
x^(K+1)/(K!*(2K+1)*N^(2K+1)), K=100. Reciprocal zeta is recomputed through
exact Bernoulli numbers and Machin's pi identity, never Arb zeta. All independent
arithmetic uses 900-bit balls. This method shares no evaluation code with riesz.

The coefficient EGF is separately checked using finite differences at 1600 bits.
The bound |c_k|<=zeta(2)<2 follows from its Mobius representation. Thus 2 times
the absolute exponential-series majorant, multiplied by x*e^-x, bounds its tail.
"""
import json, sys, pathlib, hashlib, math
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 numberdb._write import to_text
from generate import RieszValues, riesz, truncation
from check import exactness, measure, names_its_rings, agrees_with
import sage.version

HERE = pathlib.Path(__file__).resolve().parent
F = RealBallField(900)
N, K = 16, 100


def decimal_rational(text):
    text = str(text)
    mant, *ex = text.lower().split('e')
    exponent = int(ex[0]) if ex else 0
    if '.' in mant:
        whole, frac = mant.split('.')
        return QQ(int(whole+frac)) * QQ(10)**(exponent-len(frac))
    return QQ(int(mant)) * QQ(10)**exponent


def written_bounds(text):
    text = str(text)
    if '.' not in text and 'e' not in text.lower():
        q = QQ(text); return q,q,QQ(0)
    mant, *ex = text.lower().split('e')
    exponent = int(ex[0]) if ex else 0
    ulp = QQ(10)**(exponent - (len(mant.split('.')[1]) if '.' in mant else 0))
    centre = decimal_rational(text)
    return centre-ulp, centre+ulp, ulp


def enclosed(text, value):
    lo, hi, ulp = written_bounds(text)
    if ulp == 0:
        return value == lo
    return bool(value.is_finite() and F(lo) < value and value < F(hi))


def mobius(n):
    sign, p = 1, 2
    while p*p <= n:
        if n % p == 0:
            n //= p; sign = -sign
            if n % p == 0:
                return 0
        p += 1
    return -sign if n > 1 else sign


def atan_inverse(q):
    # Alternating rational terms with an explicit next-term error.
    total=F(0); j=0
    target=QQ(1)/ZZ(10)**285
    while True:
        total += F(QQ((-1)**j) / ((2*j+1)*ZZ(q)**(2*j+1)))
        j+=1
        tail=QQ(1)/((2*j+1)*ZZ(q)**(2*j+1))
        if tail < target:
            return total + F(0,tail)


def bernoullis(last):
    B=[QQ(1)]
    for n in range(1,last+1):
        B.append(-sum((QQ(math.comb(n+1,j))*B[j] for j in range(n)),QQ(0))/QQ(n+1))
    return B


PI=16*atan_inverse(5)-4*atan_inverse(239)
B=bernoullis(2*K)
MU=[0]+[mobius(n) for n in range(1,N+1)]
assert MU[1:7] == [1,-1,-1,0,-1,1]
assert B[:5] == [QQ(1),-QQ(1)/2,QQ(1)/6,QQ(0),-QQ(1)/30]
assert PI.is_finite() and PI.accuracy()>880 and PI>3 and PI<4
RECIP=[2*F(math.factorial(2*k+2))/((2*PI)**(2*k+2)*F(abs(B[2*k+2]))) for k in range(K)]
RESIDUAL=[RECIP[k]-F(sum((QQ(MU[n])/ZZ(n)**(2*k+2) for n in range(1,N+1)),QQ(0))) for k in range(K)]


def independent(x):
    x=QQ(x)
    if x==0: return ZZ(0)
    small=sum((F(QQ(MU[n])/n**2)*(-F(x)/n**2).exp() for n in range(1,N+1)),F(0))
    tailpoly=F(0)
    for k in reversed(range(K)):
        tailpoly=tailpoly*(-F(x))+RESIDUAL[k]/math.factorial(k)
    remainder=x**(K+1)/QQ(math.factorial(K)*(2*K+1)*N**(2*K+1))
    return F(x)*(small+tailpoly)+F(0,remainder)


def egf_values(xs):
    E=RealBallField(1600)
    M,_=truncation(QQ(max(xs)),135)
    diffs=[1/E(2*j+2).zeta() for j in range(M)]
    coeff=[]
    for j in range(M):
        coeff.append(diffs[0]); diffs=[diffs[k]-diffs[k+1] for k in range(len(diffs)-1)]
    results={}
    for x in xs:
        if x==0: results[x]=ZZ(0);continue
        t,s=E(1),E(0)
        for j,c in enumerate(coeff):
            s+=t*c; t*=E(x)/(j+1)
        bound=2*E(x)*t/(1-E(x)/(M+1))*(-E(x)).exp()
        value=E(x)*(-E(x)).exp()*s+E(0,bound.upper())
        assert value.is_finite() and value.accuracy()>450
        results[x]=value
    return results,coeff,M


def main():
    # Controls precede the whole-range comparison.
    assert enclosed('0.500000',F(QQ(1)/2))
    assert not enclosed('0.500000',F(QQ(500010)/10**6))
    assert (1/F(2).zeta()).overlaps(6/PI**2)
    assert not (1/F(2).zeta()).overlaps(6/PI**2+F(QQ(1)/1000))
    assert riesz(0)==0 and independent(0)==0
    assert independent(1)>0 and independent(2)<0
    generator=RieszValues(); params=list(generator.enumerate())
    assert params==[{'x':str(x)} for x in range(101)]
    values={int(p['x']):generator.value(p,100) for p in params}
    assert not exactness(values)
    assert not names_its_rings(str(HERE/'generate.py'))
    egf, coefficients, M=egf_values(list(values))
    # Independent Bernoulli inputs against the ingredients, not stored decimals.
    assert all(RECIP[k].overlaps(1/F(2*k+2).zeta()) for k in range(K))
    # Closed forms fix the signs and normalization of the EGF coefficients.
    c0=6/PI**2; c1=c0-90/PI**4; c2=c0-180/PI**4+945/PI**6
    assert all(coefficients[k].overlaps(v) for k,v in enumerate([c0,c1,c2]))
    companion=json.loads((HERE/'T644.json').read_text())['Numbers']
    for k in range(101):
        assert enclosed(companion[str(k)],F(coefficients[k]))
    if '--stored' in sys.argv:
        document=numberdb.table('T645'); numbers=document['Numbers']
        assert list(numbers)==[str(x) for x in range(101)]
    else:
        numbers={str(x):to_text(v,100) for x,v in values.items()}
    rows=[]; independent_values={}
    for x,v in values.items():
        text=numbers[str(x)]
        if isinstance(text,dict):text=text['number']
        ind=independent(x); independent_values[x]=ind
        assert enclosed(text,v) and enclosed(text,ind) and enclosed(text,egf[x]),x
        if x:
            assert v.is_finite() and v.accuracy()>350
            assert ind.is_finite() and ind.accuracy()>600
            assert v.overlaps(ind) and v.overlaps(egf[x])
            # Cislo-Wolf (33), including the full lower-order terms.
            k=F(x)
            bound=3*PI.sqrt()/16/k**(QQ(3)/2)+F(27)/2*(-F(3)).exp()/k**2+15*PI.sqrt()/32/k**(QQ(5)/2)+128*(-F(4)).exp()/k**3
            lo,hi,ulp=written_bounds(text)
            stored_ball=F((lo+hi)/2,ulp)
            clo,chi,cu=written_bounds(companion[str(x)])
            ck=F((clo+chi)/2,cu)
            assert abs(stored_ball/k-ck)<bound
            # Every decimal is sensitive to a deliberate ten-ulp shift.
            shifted=ind+10*F(ulp)
            assert not enclosed(text,shifted)
        rows.append({'x':str(x),'number':text,'primary_accuracy_bits':int(v.accuracy()) if x else None,'independent_accuracy_bits':int(ind.accuracy()) if x else None,'primary_radius':str(v.rad()) if x else '0','independent_radius':str(ind.rad()) if x else '0','terms':truncation(QQ(x),100)[0] if x else 0,'three_enclosures_inside_written_interval':True})
    # agrees_with is intended for exact equality; apply it to independently
    # rounded rational decimal centers (all 100 significant digits agree).
    assert not agrees_with({x:decimal_rational(to_text(v,100)) for x,v in values.items()},lambda x:decimal_rational(to_text(independent_values[x],100)))
    evidence={'table':'T645','stage':'build','stored_values_checked':'--stored' in sys.argv,'skill_manifest':'agents/runs/20261101T075931Z-build-skills/manifest.json','skills':json.loads((HERE/'skill-manifest.json').read_text()),'versions':{'sage':sage.version.version,'numberdb':numberdb.__version__},'entries':len(rows),'measurement_raw_balls':measure(values),'measurement_written':measure(numbers),'primary_bits_range':[numberdb.bits(100,losing=82),numberdb.bits(100,losing=280)],'minimum_primary_accuracy_bits':min(r['primary_accuracy_bits'] for r in rows if r['x']!='0'),'minimum_independent_accuracy_bits':min(r['independent_accuracy_bits'] for r in rows if r['x']!='0'),'independent_N':N,'independent_K':K,'independent_working_bits':900,'egf_terms':M,'egf_working_bits':1600,'checks':['exact origin','known rational nonzero control','deliberate-error controls','all 100 reciprocal zeta factors from exact Bernoulli numbers and Machin pi','all 101 accelerated Mobius values','all 101 EGF values with tail bound','first three coefficient closed forms','T644 coefficient intervals k=0..100','Cislo-Wolf equation (33) for k=1..100','every primary and independent enclosure contained in written last-place interval','all 100 ten-ulp perturbations rejected'], 'rows':rows}
    print(json.dumps(evidence,indent=2))


if __name__=='__main__':main()