对素幂数$p^k$（$p$是素数，$k\in\mathbb N$）有
$$\sigma_2(p^k)=1+p+p^2+\ldots+p^k.$$
$\sigma_2(n)$是积性函数，因此若$n=p_1^{a_1}\ldots p_k^{a_k}$是n的质因数分解，则
$$\sigma_2(n)=\prod_{i=1}^k\sigma_2(p_i^{a_i})=\prod_{i=1}^k\left(\sum_{j=0}^{a_i}p_i^j\right).$$
接下来没有什么太好的手段，只有逐一计算。
先用筛法生成$\sqrt{64000000}=8000$以内的所有素数，之后用同样的思想，建一个列表s用来存放计算结果，初值赋为1。
遍历素数表，对每一个素数p，遍历其所有倍数n并计算n的质因数分解中p的最大素幂因子$p^k$，并将$s(n)$更新为$s(n)\sigma_2(p^k)$。
接下来遍历所有数n，先检查$s(n)$是否为1，若$s(n)=1$，说明这是个（大于8000的）素数，遍历其倍数并执行上述操作。然后检查$s(n)$是否为平方数再求和即可。
结果是$1922364685$。
注：以下为Python 3代码，因math.isqrt是Python 3的新函数。本题中的整数非常大，无法转换为标准类型中的浮点数，因此标准库的math.sqrt无法用于计算。若使用Python 2，可以用标准库的decimal包开平方，或手写一个开平方（比如用牛顿法）的函数。本题的复杂度高，因此代码中输出了一部分运行进度以防假死。
注2：$\sigma_k$函数及其积性一般在解析数论或竞赛课程中涉及。
import math
target=8000
d=[0]*target
n=2
while n < target:
for i in range(n*n,target,n):
d[i]=d[i]+1
n=n+1
while n < target and d[n]>0:
n=n+1
primes=[k for k in range(2,target) if d[k]==0]
print("prime list generated")
bound=64000001
s=[1]*bound
for i,p in enumerate(primes):
pk=[1]
ps=[1]
while pk[-1] < bound:
pk.append(pk[-1]*p)
ps.append(ps[-1]+pk[-1]*pk[-1])
for n in range(p,bound,p):
k=1
while n%pk[k+1]==0:
k+=1
s[n]*=ps[k]
if i%100==0:
print(i,p,"checked")
r=1
for n in range(2,bound):
if s[n]==1:
q=1+n*n
for x in range(n,bound,n):
s[x]*=q
m=math.isqrt(s[n])
if m*m==s[n]:
r+=n
if n%1000000==0:
print(n,r,"completed")
print(r)
|
