首先用筛法生成一定范围内的素数，这一范围可以通过多次实验合理选择，由于撰写本文时已经算出了结果，因此这里选择上界为二十万，大约有一万八千个素数。
接下来需要有序生成平立方数，这可以用一个堆来实现，堆的长度始终与上一步素数表的长度一致，堆中元素的形式为$(z,i,j)$，其中$z=p_i^3p_j^2$，$p_i,p_j$分别是第i个和第j个素数。
每次弹出堆顶元素$(z,i,j)$，就是下一个平立方数，然后只要$i\neq j+1$，就把$(p_i^3p_{j+1}^2,i,j+1)$压入堆，否则压入$(p_i^3p_{j+2}^2,i,j+2)$。
最后需要判断是否是强合数，这涉及到数字的变化和素性检验。
注意到改变一个数的某位数字就是把这个数加上$d\cdot 10^k$，其中$d$是$[-9,9]$之间的整数，例如把231变为271就是把231加上40。
因此若z能被2或5整除，那么改变任意十位及以上的数字所得都还是合数，只需检验将个位数变为1、3、7、9的情况（事实上大多数强合数都是2或5的倍数，见OEIS A118118）。
若z不能被2或5整除，则再逐一改变其各位数字并作素性检验。
素性检验采用混合模式，由于已经生成了二十万以内的素数，因此小于二十万的数直接查表，二十万到四百亿之间的数用该表中的素数逐个作整除检验，超过四百亿的用Miller-Rabin。
最终结果是$229161792008=2^3\cdot 169249^2$。
注：出乎本人意料，本题初看计算量惊人，但实际程序几乎瞬间就给出了结果。
import heapq
target=200000
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
primes=[k for k in range(2,target) if d[k]==0]
primeSet=set(primes)
n=len(primes)
seeds=[2,325,9375,28178,450775,9780504,1795265022]
def quickPower(x,y,q):
if y <= 5:
return (x**y)%q
else:
z=quickPower(x,y/2,q)
if y%2:
return (z*z*x)%q
else:
return (z*z)%q
def millerRabin(m):
if m <= target:
return m in primeSet
elif m <= target*target:
k=0
while k < n and primes[k]*primes[k] <= m and m%primes[k]:
k+=1
return k>=n or m%primes[k]>0
elif m%2==0:
return False
else:
x=m-1
while x%2==0:
x=x>>1
for s in seeds:
a=quickPower(s,x,m)
c=x
while c < m-1 and a!=1 and a!=m-1:
a=(a*a)%m
c=c << 1
if c>m-1 or (a==1 and c>x):
return False
else:
a=(a*a)%m
c=c << 1
while c <= m-1:
if a!=1:
return False
a=(a*a)%m
c=c << 1
return True
def check(z):
if any((millerRabin(z-z%10+1),millerRabin(z-z%10+3),millerRabin(z-z%10+7),millerRabin(z-z%10+9))):
return False
else:
if z%2==0 or z%5==0:
return True
else:
u=str(z)
for i in range(1,len(u)):
w=10**i
v=z-int(u[-i-1])*w
if any(millerRabin(v+j*w) for j in range(10)):
return False
return True
h=[(72,0,1)]+[(primes[i]*primes[i]*primes[i]*4,i,0) for i in range(1,n)]
heapq.heapify(h)
r=[]
while len(r) < 200:
z,i,j=heapq.heappop(h)
j+=1
if i==j:
j+=1
heapq.heappush(h,(primes[i]*primes[i]*primes[i]*primes[j]*primes[j],i,j))
if "200" in str(z) and check(z):
r.append(z)
print r[-1]
|
