from sympy import sieve N=10000000 s1=s2=0 g=[0] for k in range(1,N+1): g.append(g[-1]+s1/k) s1+=1/k s2+=1/(k*k) mu=list(sieve.mobiusrange(1,N+1)) print("{0:.4f}".format(0.5*(N-2+s1*s1-s2)+sum(mu[r-1]*g[N//r]/(r*r) for r in range(2,N+1))))