"""Independent direct-integral checks for T681; no writes.

Run: sage -python quadrature.py, with generate.py alongside it.
Mpmath tanh-sinh quadrature uses the Mobius product, not either finite
Clausen formula. This is a heuristic numerical check, NOT the certification
of the stored digits. The generator's finite Arb formulas certify those.

The integral is split at every j/n. All logarithms of sine are evaluated
with exact small residue numerators and endpoint distances, avoiding a loss
from evaluating sin(pi*m) close to an integer m. Removable zeros are split
as well, so no quadrature interval crosses a root of a Mobius factor.
"""
import json
from math import gcd
import time
import mpmath as mp
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.arith.misc import moebius
from numberdb._write import to_text
from generate import moment


def integral(n, dps):
    with mp.workdps(dps):
        terms=[(n//int(d),int(moebius(n//d))) for d in ZZ(n).divisors()
               if moebius(n//d)]
        total=mp.mpc(0)
        for j in range((n+1)//2):
            def integrand(u):
                if not u or u==1:
                    return mp.mpc(0)
                f=mp.mpf(0)
                for q,mu in terms:
                    y=(j%q)+u
                    if y>mp.mpf(q)/2:
                        # Preserve 1-u, rather than cancelling at a large q.
                        y=(q-1-j%q)+(1-u)
                    f+=mu*mp.log(2*mp.sinpi(y/q))
                return mp.mpc(f*f,f*f*f)
            top=mp.mpf(1)/2 if j==(n-1)//2 else mp.mpf(1)
            total += mp.quad(integrand,[0,top],method='tanh-sinh')
        return total*2/n


def main():
    ns=[n for n in range(1,106,2) if moebius(n)!=0]
    # Nonzero positive and negative known controls must work first.
    mp.mp.dps=60
    generating = lambda s: mp.gamma(1+s)/mp.gamma(1+s/2)**2
    gamma_second=mp.diff(generating,0,2)
    gamma_third=mp.diff(generating,0,3)
    assert abs(gamma_second-mp.pi**2/12)<mp.mpf('1e-55')
    assert abs(gamma_third+mp.mpf(3)*mp.zeta(3)/2)<mp.mpf('1e-55')
    ctrl=integral(1,50)
    assert abs(ctrl.real-mp.pi**2/12)<mp.mpf('1e-42')
    assert abs(ctrl.imag+mp.mpf(3)*mp.zeta(3)/2)<mp.mpf('1e-42')
    print('CONTROL: direct n=1 integral matches pi^2/12 and -3*zeta(3)/2 to 42 decimal places',flush=True)
    rows=[]
    for n in ns:
        started=time.monotonic()
        val=ctrl if n==1 else integral(n,50)
        targets=[mp.mpf(to_text(moment(n,k),digits=100)) for k in (2,3)]
        errors=[abs(val.real-targets[0]),abs(val.imag-targets[1])]
        assert max(errors)<mp.mpf('1e-40'), (n,errors)
        row={'n':n,'m2':mp.nstr(val.real,45),'m3':mp.nstr(val.imag,45),
             'absolute_errors':[mp.nstr(e,8) for e in errors],
             'seconds':round(time.monotonic()-started,3)}
        rows.append(row); print('QUAD',json.dumps(row),flush=True)
    print('QUADRATURE_JSON '+json.dumps({'table':'T681','dps':50,'absolute_tolerance':'1e-40',
        'method':'mpmath tanh-sinh; Mobius log-sine sum; split at j/n; combined second and third moments',
        'mpmath_version':mp.__version__,'values':rows,
        'skill_manifest':'agents/runs/20261102T162707Z-build-skills/manifest.json'}),flush=True)


if __name__=='__main__': main()