先取$d=1$,则$1+n$是素数。
再考虑n的素因子p,若$p^2|n$,则$p+n/p=p(1+n/p^2)$不是素数,从而n是无平方因子数。
因此先用筛法找出不超过$10^8$的所有素数,对每个素数p,考察$m=p-1$是否是无平方因子数,最后枚举m的约数检验是否符合条件即可。
最终结果是$1739023853137$。
注:因用到Walrus算符和math.isqrt函数,以下代码为Python 3。此外,代码中打印了进度信息。
import math
target=10**8
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
primeList=[k for k in range(2,target) if d[k]==0]
primeSet=set(primeList)
print("prime list generated with length", len(primeList))
tick=len(primeList)//100
r=0
for k,n in enumerate(primeList):
m=n-1
i=0
while (q:=primeList[i]*primeList[i])<=m and m%q>0:
i+=1
if q>m and all((m//i+i) in primeSet for i in range(2,math.isqrt(m)+1) if m%i==0):
r+=m
if k%tick==0:
print(k//tick,"percent completed")
print(r)
|