记$g(n)=\gcd(n,n')$，$g(n)$的值可参考OEIS A085731。特别地，若$n=p_1^{a_1}\ldots p_k^{a_k}$，则
$$n'=a_1\ldots a_k\cdot p_1^{a_1-1}\ldots p_k^{a_k}$$
若令
$$f(p^a)=g(p^a)-g(p^{a-1}),$$
则$f$是积性函数，从而若令$P$为所有PN数(即powerful number，与无平方因子数相反，若$n$是PN数则对其每个素因子$p$而言$p^2$都能整除$n$)，则
$$\sum_{k=2}^Ng(k)=-1+\sum_{n\in P}\left\lfloor\frac{N}{n}\right\rfloor f(n),$$
此处的求和是PN筛的标准公式（见此处）。
直接计算即得最终结果$8907904768686152599$。
注：因使用sympy库以及math.gcd等函数，以下代码为Python 3。代码中未打印进度，但运行需要数分钟。
import sympy,math
N=5*(10**15)
preLists=[]
def f(a,m,x):
res=m*x
for p in range(a, len(preLists)):
if preLists[p][0][0] > m:
break
for q,r in preLists[p]:
if q>m:
break
res+=f(p+1,m//q,x*r)
return res
for p in sympy.primerange(math.isqrt(N)):
q,r,e= p*p,1,2
pre=[]
while q<=N:
d=math.gcd(e*(q//p), q)
if d-r:
pre.append((q,d-r))
q*=p
e+=1
r=d
preLists.append(pre)
print(f(0,N,1)-1)
|
