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.rings.complex_arb import ComplexBallField
from sage.arith.misc import moebius

def riemann_sum_error(N, digits=100):
    N=ZZ(N)
    if N<1: raise ValueError('N must be positive')
    if N<=2: return ZZ(2)
    R=RealBallField(numberdb.bits(digits,losing=120))
    C=ComplexBallField(R.precision())
    M=[sum((moebius(j) for j in range(1,n+1)),ZZ(0)) for n 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())==0
    jumps={QQ(0):QQ(0),QQ(1):QQ(0)}
    for k,v in w.items():
        for j in range(1,k):
            a=QQ(j)/k; jumps[a]=jumps.get(a,QQ(0))-v
    cuts=sorted(jumps); c=QQ(0); norm=R(0)
    for a,b in zip(cuts,cuts[1:]):
        c+=jumps[a]
        if c: norm+=c*c*(C(b).psi()-C(a).psi()).real()/N
    cross=R(0)
    for k,v in w.items():
        L=R(QQ(k)/N).log()
        cross+=R(v*QQ(k)/N)*((1-R.euler_constant())*L-L*L/2)
    return 2+norm-2*cross

# Stored coverage: every N=3,...,100; the function accepts all positive N.
riemann_sum_error(100)