若$\varphi(n)=13!$，且p是n的一个素因子，则$p-1$能整除$13!$。
$$13!=6227020800=2^10\times 3^5\times 5^2\times 7\times 11\times 13,$$
因此先找出其所有因数d，共有$11\times 6\times 3\times 2\times 2\times 2$大约1600个，并从中找出满足$d+1$是素数的因子，这样的素数共有459个。
接下来作深度优先搜索，在栈中保存三元组当前数，当前数的Totient函数值，当前数的最大素因子下标。
每次弹出栈顶元素后从其最大素因子（因为存在素幂因子）开始遍历素数表并计算相应的Totient函数值，把值等于$13!$放入结果列表，值小于且能整除$13!$的压入栈。
最后排序结果列表可得答案是$23507044290$。
注：计算量较大，代码中打印了进度信息。
target=80000
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)
def isPrime(x):
if x <= target:
return x in primeSet
else:
i=0
while i < len(primeList) and x%primeList[i]>0:
i+=1
return i==len(primeList)
primeCandidates=[]
for a2 in range(11):
for a3 in range(6):
for a5 in range(3):
for a7 in range(2):
for a11 in range(2):
for a13 in range(2):
d=(2**a2)*(3**a3)*(5**a5)*(7**a7)*(11**a11)*(13**a13)
if isPrime(d+1):
primeCandidates.append(d+1)
primeCandidates.sort()
print "primes prepared with length", len(primeCandidates)
M=6227020800
r=[]
q=[(p,p-1,i) for i,p in enumerate(primeCandidates)]
while q:
n,x,k=q.pop()
i=k
while i < len(primeCandidates):
p=primeCandidates[i]
y=x*p if n%p==0 else x*(p-1)
if y>M:
break
elif y==M:
r.append(n*p)
if len(r)%10000==0: print len(r), "collected"
elif M%y==0:
q.append((n*p,y,i))
i+=1
r.sort()
print r[149999]
|
