"""T711: Baez-Duarte R_N squared logarithmic-target errors.

Install with `sage -pip install numberdb`.
Run `sage -python generate.py` to verify the
entire stored table, or `sage -python generate.py --publish` to fill its draft.
The latter writes entries only, never offers or publishes the table itself.

N is an integer >=1. Default coverage is every 3<=N<=100. N=1,2 are zero
approximants with exact error 2, retained as controls rather than generic rows.
The cutoff is a computational extent, not part of the mathematical domain.
Uses exact Mobius/Mertens weights and the finite rational autocorrelation
formula (BBLS, arXiv:math/0306251, Prop.89). The logarithmic cross term follows
by Mellin integration; its common Stieltjes term cancels by sum M(N/k)=0,
k=1,...,N-1. See derivation.md and independent_check.py.
All transcendental evaluations use Arb balls at bits(digits)+120 bits.
Real logarithms and positive pi are used; there is no infinite numerical sum.
The run evidence records the measured loss and independent comparisons.
"""
import argparse
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.arith.misc import moebius

@lru_cache(None)
def mertens(n):
    return sum((ZZ(moebius(j)) for j in range(1,n+1)),ZZ(0))

@lru_cache(None)
def cotangents(q,bits):
    R=RealBallField(bits)
    return tuple((R.pi()*R(QQ(m)/q)).cot() for m in range(1,q))

@lru_cache(None)
def vasyunin(p,q,bits):
    R=RealBallField(bits)
    return sum((R(QQ((m*p)%q)/q)*c for m,c in enumerate(cotangents(q,bits),1)),R(0))

@lru_cache(None)
def autocorrelation(r,bits):
    R=RealBallField(bits); p=int(r.numerator());q=int(r.denominator()); x=R(r)
    return (1-x)*x.log()/2+(1+x)*((2*R.pi()).log()-R.euler_constant())/2-R.pi()*(vasyunin(p,q,bits)+vasyunin(q,p,bits))/(2*q)

def components(n,digits=100):
    n=ZZ(n)
    if n<1:raise ValueError('N must be >= 1')
    bits=numberdb.bits(digits,losing=120); R=RealBallField(bits)
    terms=[(k,mertens(int(n//k))) for k in range(1,int(n)) if mertens(int(n//k))]
    assert sum((m for k,m in terms),ZZ(0))==0
    weights=[(k,QQ(m)/k) for k,m in terms]
    norm=R(0)
    for i,(j,wj) in enumerate(weights):
        for k,wk in weights[:i+1]:
            norm+=(1 if j==k else 2)*R(wj*wk*QQ(j)/n)*autocorrelation(QQ(k)/j,bits)
    cross=R(0)
    for k,m in terms:
        logtheta=R(QQ(k)/n).log()
        cross+=R(QQ(m)/n)*((1-R.euler_constant())*logtheta-logtheta**2/2)
    return norm,cross

class RiemannErrors(numberdb.Generator):
    table='T711'
    parameters=('N',)
    type='R'
    digits=100
    rigour='proven'
    def enumerate(self):
        for n in range(3,101):yield {'N':n}
    def value(self,params,digits):
        n=ZZ(params['N'])
        if n<1:raise ValueError('N must be >= 1')
        if n<=2:return ZZ(2)
        norm,cross=components(n,digits)
        result=2+norm-2*cross
        assert result.is_finite() and result>0
        assert result.accuracy()>=numberdb.bits(digits)+30
        return result

if __name__=='__main__':
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--publish',action='store_true')
    args=parser.parse_args()
    gen=RiemannErrors()
    if args.publish:print(gen.publish())
    else:print(gen.verify(sample=None))