根据OEIS A034676提供的信息，$S(n)$是积性函数，而显然$S(p^e)=p^{2e}+1$，因此先用筛法找出不超过$10^8$的素数，再找出每个素数在$10!$的质因数分解中的幂，汇总即得结果$98792821$。
注：因数据较大，为了更好使用内存和利用内建函数作快速幂，以下代码为Python 3。
target=10**8
mod=10**9+9
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
print("prime list generated, start computing...")
res=1
for p in range(2,target):
if d[p]==0:
q=p
e=0
while q <= target:
e+=target//q
q*=p
res=(res*(pow(p,2*e,mod)+1))%mod
print(res)
|
