先考虑n模不超过10的整数的情况：
2到4：$n-3$是素数，因此n是偶数，且3不能整除n。从而4需要能整除n，否则$n\pm 4$和$n\pm 8$的因子中无4，也就不能表示4，而不能成为实用数。
5到7：$n\pm 3, n\pm9$都是素数，因此n需要是5的倍数，同理n模7只能余0,1,6之一。模6的情况则已经蕴含在模2和模3之中。
8：由上文分析可知1,2,4都是n的约数，从而也是$n\pm 4$和$n\pm 8$的约数。5是n的约数，因此不是$n\pm 4$和$n\pm 8$的约数。若8能整除n，则其不能整除$n\pm 4$，而这两者中至少有一个无法被3整除，因而此时$n\pm 4$中至少有一个的约数和无法表示8，因此n模8只能余4。
综上，由中国剩余定理可得n模$840(=3\times 5\times 7\times 8)$余$\pm 20$。
枚举这样的n并检验，直到获得四个工程师天堂数为止。
作素性检验需要的素数上限难以预先判断，经测试，生成五万以内的素数已足够。此外，需要注意此处要求的是连续是素数，因此不仅需要检验$n\pm 9$和$n\pm 3$的素性，还需要验证从$n-9$到$n+9$之间的其他数不是素数。
对于实用数的判断，注意到若$m=p_1^{a_1}\ldots p_k^{a_k}$是实用数$n$的质因数分解，则不难归纳出对$i=1,\ldots,k$其都需要满足（或参考此处）:
$$p_i\le 1+\prod_{j=1}^{i-1} \dfrac{p_j^{e_j+1}-1}{p_j-1}。$$
最终结果是$2039506520$。对应的四个数是219869980, 312501820, 360613700, 1146521020。
注：计算量较大，代码中打印了进度信息。
target=50000
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)
def isPractical(n):
m=n
i=0
r=1
while i < len(primeList) and m>1:
if m%primeList[i]==0:
if primeList[i]>r+1:
return False
else:
p=primeList[i]
while m%primeList[i]==0:
m/=primeList[i]
p*=primeList[i]
r*=(p-1)/(primeList[i]-1)
i+=1
return m==1 or m <= r+1
paradise=[]
n=820
while len(paradise) < 4:
if isPrime(n-9) and isPrime(n-3) and isPrime(n+3) and isPrime(n+9) and not any(isPrime(n+k) for k in [-8,-7,-6,-5,-4,-2,-1,0,1,2,4,5,6,7,8]) and isPractical(n-8) and isPractical(n-4) and isPractical(n) and isPractical(n+4) and isPractical(n+8):
paradise.append(n)
print paradise
n+=40
if isPrime(n-9) and isPrime(n-3) and isPrime(n+3) and isPrime(n+9) and not any(isPrime(n+k) for k in [-8,-7,-6,-5,-4,-2,-1,0,1,2,4,5,6,7,8]) and isPractical(n-8) and isPractical(n-4) and isPractical(n) and isPractical(n+4) and isPractical(n+8):
paradise.append(n)
print paradise
n+=800
if (n+20)%8400000==0:print (n+20)/8400000,"w pairs checked"
print sum(paradise)
|
