本题数据量很大，直接用筛法生成所有数的质因数分解在时间和空间上都不可取，因此先只用筛法生成四千万以内的所有素数，同时放在一个列表（用于循环）和一个集合（用于查询）中。
用一个列表缓存中间计算所得的数的数链长度，依次遍历这些素数，每个数只迭代最多24次，并记录中间所得的数。
每次迭代先查表，如表中有结果直接和本轮迭代所得的数列合并记录。
如表中无结果，就用公式
$$\varphi(n)=n\prod_p\frac{p-1}{p},$$
计算欧拉totient函数的值，此处$p$取遍$n$的所有素因子。
具体过程是先令临时变量$m=n$，以计算n的质因数分解，n则可以直接逐步更新为$\varphi(n)$。
遍历素数表，若某个素数p能整除m，就把n更新为$n(p-1)/p$，并不断把m更新为$m/p$直到p不再是m的约数。
当$m=1$或为素数（所以在第一步中需要维护素数集）时本次迭代终止，若m为素数，n还需更新为$n(m-1)/m$。
最终结果是$1677366278943$。
注：本题迭代中的记忆化搜索和素数集查询都不可省略，否则将极大增加时间复杂度。为防假死，以下代码中还输出了运行进度信息。
target=40000001
bound=6325
d=[0]*target
n=2
while n < bound:
for i in range(n*n,target,n):
d[i]+=1
n+=1
while n < bound and d[n]>0:
n+=1
primes=[k for k in range(2,target) if d[k]==0]
primeSet=set(primes)
print "prime set generated"
chain=[0]*target
chain[1]=1
r=0
for j,p in enumerate(primes):
d=[p]
n=p-1
while chain[n]==0 and len(d) < 25:
d.append(n)
m=n
i=0
while m>1 and m not in primeSet:
if m%primes[i]==0:
n=(n/primes[i])*(primes[i]-1)
while m%primes[i]==0:
m/=primes[i]
i+=1
if m>1:
n=(n/m)*(m-1)
if chain[n]>0:
for i,q in enumerate(d):
chain[q]=min(chain[n]+len(d)-i,26)
else:
for q in d:
chain[q]=26
if chain[p]==25:
r+=p
if j%10000==0:
print "Progress:",j/10000,"/244; Current r:", r
print r
|
