$n^2$的最后一位数是周期序列$0,1,4,9,6,5,6,9,4,1$。因此：
若$n^2+1$是素数，那么$n^2$的最后一位数字不能是奇数（否则$n^2+1$将是偶数）或4（否则$n^2+1$将能被5整除）。
若$n^2+9$是素数，那么$n^2$的最后一位数字不能是6，否则$n^2+9$将能被5整除。
因此$n^2$的最后一位数只能是0，即n是10的倍数。
现考虑$n \bmod 7$：
$n \bmod 7\neq0$，否则$n^2+7$将能被7整除而不是素数。
$n \bmod 7\neq1$，否则$n^2+13$将能被7整除而不是素数。
$n \bmod 7\neq2$，否则$n^2+3$将能被7整除而不是素数。
$n \bmod 7\neq5$，否则$n^2+3$将能被7整除而不是素数。
$n \bmod 7\neq2$，否则$n^2+13$将能被7整除而不是素数。
因此$n \bmod 7$只能为3或4。类似地，n不能被3整除，否则$n^2+3$将能被3整除而不是素数。
综上所述，由中国剩余定理可得$n=10,80,130,200 ,\pmod{210}$，从而将需要检验的n减少到了三百万个左右。
此外，题中要求连续素数，因此还需检验$n^2+5$, $n^2+11$, $n^2+15$, $n^2+17$, $n^2+19$, $n^2+21$, $n^2+23$, $n^2+25$，这些数应为合数。
最后，由于$n^2$可以大至$10^{16}$，单纯的筛法无论在空间和时间上都不太合适，因而此处混合使用素性检验，先用筛法生成一百万以下的素数表，$10^6$以下直接查表，$10^6$到$10^{12}$之间的用该表中的素数依次检验整除性，$10^{12}$以上的使用Miller-Rabin法。
结果是$676333270$。
seeds=[2,325,9375,28178,450775,9780504,1795265022]
p1000=[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997]
smallPrimes=set(p1000)
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 <=1000:
return m in smallPrimes
elif m <=1000000:
k=0
while k < len(p1000) and p1000[k]*p1000[k] <=m and m%p1000[k]:
k+=1
return k>=len(p1000) or m%p1000[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
r=0
for i in range(714287):
for j in (10,80,130,200):
n=210*i+j
if n < 150000000:
z=n*n
if all(millerRabin(z+k) for k in (1,3,7,9,13,27)) and all(not millerRabin(z+k) for k in (5,11,15,17,19,21,23,25)):
r+=n
if i%10000==0 and j==10: print i,n,r
print r
|
