"""Dual Mobius squared approximation errors -- numberdb.org/T709.

sage -pip install numberdb
sage -python generate.py              # verify only
sage -python generate.py --publish    # fill authorized target, never offer
For a private draft configure NUMBERDB_API_KEY. N is a positive integer;
the default enumeration is every N=1,...,30. See derivation.md for the
source translation, the finite cotangent formula and the enclosure argument.
All arithmetic weights are exact. Arb enclosures use 128 guard bits beyond
the requested digits; the computed prefix has separately measured accuracy.
N=1 has exact error 1 because the dual approximant vanishes.
The gram_system component is reused from T706/T641, original T641 SHA-256
7e155af1d4a4963f874e0c536051c1cf0b7afd01e21b93d64a674c7fab8592c8.
No sibling-table decimal is an input. The direct-residual checker avoids
this Gram implementation and independently constructs the arithmetic weights.
"""
import os
import sys
from functools import lru_cache
from math import gcd
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.arith.misc import moebius

MAX_N=30
WORKING_GUARD=128

def gram_system(limit, precision):
    """Fresh finite cotangent sums, Gram entries, and target inner products."""
    R = RealBallField(precision)
    C = (2*R.pi()).log() - R.euler_constant()

    @lru_cache(None)
    def vasyunin(p, q):
        return sum((R(QQ((m*p) % q)/QQ(q)) *
                    (R.pi()*R(QQ(m)/QQ(q))).cot()
                    for m in range(1, q)), R(0))

    @lru_cache(None)
    def autocorrelation(p, q):
        x = R(QQ(p)/QQ(q))
        return (((1-x)*x.log() + (1+x)*C)/2
                - R.pi()*(vasyunin(p,q)+vasyunin(q,p))/(2*q))

    G = [[R(0) for k in range(limit)] for j in range(limit)]
    for j in range(1, limit+1):
        for k in range(1, j+1):
            g = gcd(j,k)
            G[j-1][k-1] = G[k-1][j-1] = autocorrelation(j//g,k//g)/j
    b = [(R(n).log()+1-R.euler_constant())/n for n in range(1,limit+1)]
    return R, G, b


def coefficients(N):
    """M(N/k)-M(N/(k+1))-delta(k,1), from equation (1.19)."""
    M=[ZZ(0)]
    for n in range(1,N+1):M.append(M[-1]+moebius(n))
    w=[M[N//k]-M[N//(k+1)] for k in range(1,N+1)]
    w[0]-=1
    assert sum((k*w[k-1] for k in range(1,N+1)),ZZ(0))==0
    return w

@lru_cache(None)
def computation(limit,digits):
    R,G,_=gram_system(limit,numberdb.bits(digits,losing=WORKING_GUARD))
    values=[ZZ(1)]
    for N in range(2,limit+1):
        w=coefficients(N)
        b=[R(QQ(k)/QQ(N))*(R(QQ(N)/QQ(k)).log()+1-R.euler_constant())
           for k in range(1,N+1)]
        # Dual Gram entry is j*k/N times the integer-dilation Gram entry.
        v=R(1)+2*sum((w[k-1]*b[k-1] for k in range(1,N+1)),R(0))
        v+=sum((w[j-1]*w[k-1]*R(QQ(j*k)/QQ(N))*G[j-1][k-1]
                for j in range(1,N+1) for k in range(1,N+1)),R(0))
        assert v.is_finite() and v>0
        assert v.accuracy()>=numberdb.bits(digits)+16
        values.append(v)
    return R,G,values

class DualErrors(numberdb.Generator):
    table='T709'
    parameters=('N',)
    type='R'
    digits=100
    rigour='proven'
    files=('generate.py','checks.py','integral_kernel.py','argument-selection.md',
           'derivation.md','numerical-evidence.json','skill-manifest.json','research.md',
           'verify_stored.py')

    def enumerate(self,limit=MAX_N):
        for N in range(1,limit+1):yield {'N':str(N)}

    def value(self,params,digits):
        N=int(params['N'])
        if N<1:raise ValueError('N must be a positive integer')
        if N==1:return ZZ(1)
        return computation(max(MAX_N,N),digits)[2][N-1]

if __name__=='__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        numberdb.configure(api_key=sys.stdin.read().strip())
    generator=DualErrors()
    if '--publish' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1':
        print(generator.publish(message=os.environ.get('NUMBERDB_EDIT_MESSAGE','Recompute dual Mobius errors for N=1,...,30')))
    else:
        report=generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)