设$n=p_1^{a_1}\ldots p_k^{a_k}$是阿喀琉斯数n的质因数分解，其中$p_1 < \ldots < p_k$都是素数。
显然$a_1,\ldots,a_k\ge 2$，否则n将有无平方素因子，$\gcd(a_1,\ldots,a_k)=1$，否则n将是另一数的幂。
若n是强阿喀琉斯数，那么还应有$a_k\ge 3$，否则$p_k$将是$\varphi(n)$的无平方素因子。由此可得$p_k$不超过$\sqrt[3]{10^{18}/4}$约63万。除此以外似乎也没有其他可以直接利用的性质。
先用筛法找出不超过63万的素数，对每个素数p，计算出并缓存p-1的质因数分解。
接下来从大到小依次枚举n的最大素因子，用深度优先搜索的方式对每个素因子，添加所有可能的幂并检验，最终结果是$1170060$。
注：因需要使用math.gcd和math.prod函数，以下为Python 3代码，且代码中打印了进度。
import math,functools
bound=10**18
target=629961
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]
primeDivisors={}
for k,p in enumerate(primeList):
q=p-1
r={}
i=0
while q>1:
if q%primeList[i]==0:
j=0
while q%primeList[i]==0:
q//=primeList[i]
j+=1
r[i]=j
i+=1
primeDivisors[k]=r
print("initialization completed")
def getVal(d):
return math.prod(primeList[j]**d[j] for j in d)
def isAchilles(d):
return len(d)>=2 and functools.reduce(math.gcd,d.values())==1 and all(v>=2 for v in d.values())
def getTotient(d):
t={k:d[k]-1 for k in d if d[k]>1}
for k in d:
for x in primeDivisors[k]:
t[x]=t.get(x,0)+primeDivisors[k][x]
return t
r=0
for i in range(1,len(primeList)):
k=3
p=primeList[i]**3
q=[]
while 4*p < bound:
q.append({i:k})
p*=primeList[i]
k+=1
while q:
d=q.pop()
if isAchilles(d) and isAchilles(getTotient(d)):
r+=1
if r%10000==0:print(r//10000,"W collected")
for j in range(min(d.keys())):
m=getVal(d)*primeList[j]*primeList[j]
k=2
while m < bound:
d2=d.copy()
d2[j]=d2.get(j,0)+k
q.append(d2)
m*=primeList[j]
k+=1
print(r)
|
