设b是n的素因子中模3余1型素数的数量，根据OEIS A060839中的公式，若n是9的倍数，则$C(n)+1=3^{b+1}$，否则$C(n)+1=3^{b}$。
因此本题中b需要等于4或5，最小的三个模3余1型的素书是7,13,19，因此先用筛法找出不超过$10^{11}/(9\cdot7\cdot 13\cdot19)\approx6426322$且模3余1的素数，放在一个列表中，为避免下标越界，以下代码中实际使用的上界是6500000。
用这些素数再作一次筛法，找出不超过6500000且不能被这些素数整除的所有数，放在集合q中。
接下来按照n是否为9的倍数分别枚举这些素数中四个或五个数的乘积，对每个乘积，枚举其倍数并检验倍数的因子是否只在q中，最后汇总得结果是$8495585919506151122$。
注：计算量较大，代码中打印了进度信息。
target=6500000
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 and k%3==1]
print "prime list generated"
d=[0]*target
for p in primeList:
for i in range(p,target,p):
d[i]=1
q=set([k for k in range(1,target) if d[k]==0])
print "set q generated"
bound=10**11
r=0
a=0
while primeList[a]*primeList[a+1]*primeList[a+2]*primeList[a+3]*primeList[a+4]<=bound:
b=a+1
while primeList[a]*primeList[b]*primeList[b+1]*primeList[b+2]*primeList[b+3]<=bound:
c=b+1
while primeList[a]*primeList[b]*primeList[c]*primeList[c+1]*primeList[c+2]<=bound:
d=c+1
while primeList[a]*primeList[b]*primeList[c]*primeList[d]*primeList[d+1]<=bound:
e=d+1
while primeList[a]*primeList[b]*primeList[c]*primeList[d]*primeList[e]<=bound:
m=primeList[a]*primeList[b]*primeList[c]*primeList[d]*primeList[e]
k=1
while m*k<=bound:
j=k
if j%9!=0:
for x in (primeList[a],primeList[b],primeList[c],primeList[d],primeList[e]):
while j%x==0:
j/=x
if j in q:
r+=m*k
k+=1
e+=1
d+=1
c+=1
b+=1
a+=1
print "progress: a=",a
print "Case 1 completed"
a=0
while 9*primeList[a]*primeList[a+1]*primeList[a+2]*primeList[a+3]<=bound:
b=a+1
while 9*primeList[a]*primeList[b]*primeList[b+1]*primeList[b+2]<=bound:
c=b+1
while 9*primeList[a]*primeList[b]*primeList[c]*primeList[c+1]<=bound:
d=c+1
while 9*primeList[a]*primeList[b]*primeList[c]*primeList[d]<=bound:
m=9*primeList[a]*primeList[b]*primeList[c]*primeList[d]
k=1
while m*k<=bound:
j=k
for x in (primeList[a],primeList[b],primeList[c],primeList[d]):
while j%x==0:
j/=x
if j in q:
r+=m*k
k+=1
d+=1
c+=1
b+=1
a+=1
print "progress: a=",a
print "Case 2 completed"
print r
|
