"""xexact.py -- exact rational values of e4(D), I(D), X(D) over the full period lcm(1..D).
Verifies the exact values quoted in Section 5: X(7)=X(6)=-1/324,
X(10)=...=X(13)=-217/50625, Delta X(14)=-3/4802, X(15)=822653/81033750.
PhiTilde_D is piecewise linear on (n,n+1) with slope S = sum mu(d)/d, so all
moments are finite rational sums; sigma_q is piecewise constant on (j/q,(j+1)/q).
"""
from fractions import Fraction as F
from math import lcm, gcd

def mu_sieve(n):
    mu=[0]*(n+1); mu[1]=1
    primes=[]
    is_c=[False]*(n+1)
    for i in range(2,n+1):
        if not is_c[i]: primes.append(i); mu[i]=-1
        for p in primes:
            if i*p>n: break
            is_c[i*p]=True
            if i%p==0: mu[i*p]=0; break
            else: mu[i*p]=-mu[i]
    return mu

def exact_e4(D, MU):
    L=lcm(*range(1,D+1))
    S=sum(F(MU[d],d) for d in range(1,D+1))
    half=F(1,2)
    # A_n = sum mu(d)*({n/d}-1/2)
    m2=F(0); m4=F(0)
    for n in range(L):
        # A = sum_d mu(d)*({n/d} - 1/2); on (n,n+1), PhiTilde = A + t*S
        A=sum(F(MU[d]*(n%d),d) for d in range(1,D+1)) - half*sum(MU[d] for d in range(1,D+1) if MU[d])
        # exact integrals of (A+tS)^2 and (A+tS)^4 over t in [0,1]
        m2 += A*A + A*S + S*S/3
        A2=A*A; S2=S*S
        m4 += A2*A2 + 2*A2*A*S + 2*A2*S2 + A*S*S2 + S2*S2/5
    m2/=L; m4/=L
    return m4 - 3*m2*m2, m2

def eps4(q, MU):
    divs=[t for t in range(1,q+1) if q%t==0 and MU[t]!=0]
    s2=F(0); s4=F(0)
    for j in range(q):
        v=F(0)
        for t in divs:
            # B1 at t*(j+1/2)/q  => value = frac - 1/2, frac of (t*(2j+1))/(2q)
            num=t*(2*j+1); den=2*q
            v += F(MU[t],t)*(F(num%den,den)-F(1,2))
        s2+=v*v; s4+=v**4
    s2/=q; s4/=q
    return s4-3*s2*s2

def mq(q,Y,MU):
    # sum_{m<=Y,(m,q)=1} mu(m)/m
    return sum(F(MU[m],m) for m in range(1,int(Y)+1) if gcd(m,q)==1)

def X(D, MU):
    e4,m2 = exact_e4(D,MU)
    I=F(0)
    for q in range(1,D+1):
        if q>1 and MU[q]==0: continue
        if q==1:
            I += mq(1,D,MU)**4 * eps4_1
        else:
            I += mq(q,D//q,MU)**4 * eps4(q,MU)
    return e4-I, e4

MU=mu_sieve(400)
eps4_1=F(-1,120)
for D in [6,7,10,11,12,13,14,15]:
    x,e4=X(D,MU)
    print(D, x, float(x), ' e4=',float(e4))