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

def dual_error(N, digits=100):
    N=ZZ(N)
    if N<1:raise ValueError('N must be positive')
    if N==1:return ZZ(1)
    R=RealBallField(numberdb.bits(digits,losing=128))
    w=[ZZ(0)]*N
    for m in range(1,N+1):w[N//m-1]+=moebius(m)
    w[0]-=1
    cuts=sorted({QQ(a)/k for k in range(1,N+1) for a in range(k+1)})
    total=R(0)
    for left,right in zip(cuts,cuts[1:]):
        mid=(left+right)/2
        c=1-sum(w[k-1]*(k*mid).floor() for k in range(1,N+1))
        total+=c*c*(R(1+right).psi()-R(1+left).psi())
        start=max(left,QQ(1)/N)
        if right>start:total+=c*c*(1/R(start)-1/R(right))
    return total/N

# Every N=1,...,30 is stored; this function also accepts larger positive N.
dual_error(30)