0058 螺旋素数
* * * *
拉格朗日计划
* * * *
螺旋素数

从1开始按照逆时针方向摆放自然数,可以构造出如下边长为7的螺旋数阵: $$\begin{matrix} \textcolor{red}{37} & 36 & 35 & 34 & 33 & 32 & \textcolor{red}{31} \\ 38 & \textcolor{red}{17} & 16 & 15 & 14 & \textcolor{red}{13} & 30 \\ 39 & 18 & \textcolor{red}{5} & 4 & \textcolor{red}{3} & 12 & 29 \\ 40 & 19 & 6 & 1 & 2 & 11 & 28 \\ 41 & 20 & \textcolor{red}{7} & 8 & 9 & 10 & 27 \\ 42 & 21 & 22 & 23 & 24 & 25 & 26 \\ \textcolor{red}{43} & 44 & 45 & 46 & 47 & 48 & 49 \end{matrix}$$ 对角线上共有13个数,其中8个是素数(已标记为红色),比例达到$8/13\approx62\%$。

在这个方阵外再增加一层,可得一边长为9的螺旋方阵。不断重复这一过程,当对角线上素数的比例第一次低于10%时,该螺旋数阵的边长是多少?

本题难度:



解答

边长为n时,四个角落的数分别是$n^2$、$n^2-(n-1)$、$n^2-2(n-1)$、$n^2-3(n-1)$,因此只需不断测试和统计后三项是否为素数即可,结果是$26241$。

target=100000
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)

def isPrime(m):
    if m <=100000:
        return m in primeSet
    else:
        i=0
        while primeList[i]*primeList[i] <= m and m%primeList[i]:
            i+=1
        return m%primeList[i]>0

c=8
d=13
n=7
while 10*c>=d:
    n+=2
    d+=4
    c+=(isPrime(n*n-n+1)+isPrime(n*n-2*n+2)+isPrime(n*n-3*n+3))

print n