优化埃拉斯托斯特尼筛法（查找超过一百万个素数，没有内存错误）

Question

我的代码只能处理大约一百万的范围，而我需要它处理大约一亿到十亿。

当我输入 1 到 10 亿的范围时，我不断出现内存错误。

有什么方法可以优化代码吗？

import time
from math import isqrt


def sieveprmcheck(n):
    if n < 2:
        return []
    prime_num = [True] * (n + 1)
    prime_num[0] = False
    prime_num[1] = False

    for i in range(2, isqrt(n) + 1):
        if prime_num[i]:
            for k in range(i * i, n + 1, i):
                prime_num[k] = False

    primes = [i for i in range(n + 1) if prime_num[i]]
    return primes


def segsieve():
    primes = sieveprmcheck(n)

    prime = [True] * (n - m + 1)

    for i in primes:
        lower = (m // i)

        if lower <= 1:
            lower = i + i
        elif (m % i) != 0:
            lower = (lower * i) + i
        else:
            lower = lower * i
        for j in range(lower, n + 1, i):
            prime[j - m] = False

    prime_numbers = []
    for k in range(m, n + 1):
        if prime[k - m]:
            prime_numbers.append(k)

    return prime_numbers


def palcheck(primes):
    special_numbers = set()
    for num in primes:
        if str(num) == str(num)[::-1]:  # Check if the number is a palindrome
            special_numbers.add(num)

    special_numbers_list = sorted(list(special_numbers))
    print('Number of Special Numbers:', len(special_numbers_list), ':',
          special_numbers_list[:3] + special_numbers_list[-3:])
    return special_numbers_list


if __name__ == '__main__':
    m = int(input('Please enter starting number: '))
    n = int(input('Please enter ending number: '))
    start = time.time()
    primes = segsieve()
    Special_Num = palcheck(primes)
    end = time.time()

    time_spent = end - start
    print("\nTime spent:", time_spent, "seconds")

所以它最多可以计算数亿，但我预计会从数十亿到数万亿个素数。

Answer 1

正如评论中已经建议的，您可以使用轮分解来改进筛子。

这是我不久前创建的一个示例，使用尺寸为 30 的轮子来满足您要求的使用范围。 2、3 和 5 的所有倍数均不会存储以节省内存：

import time
import numpy as np
def sieve(n):
    baseW=30
    baseP=np.array([-23,-19,-17,-13,-11,-7,-1,1])
    C=np.ones((len(baseP),len(baseP)), dtype=int)
    C1=np.zeros((len(baseP),len(baseP)), dtype=int)
    C=np.multiply(np.dot(np.diag(baseP),C),np.dot(C,np.diag(baseP)))
    C=C%baseW
    C[C>1] = C[C>1] -baseW
    for i in range(len(baseP))    :
        C1[i,:]=baseP[np.argwhere(C==baseP[i])[:,1]]
    primeB=np.ones((len(baseP),n//baseW+1), dtype=bool)
    for j in range(1,int((n**0.5-baseP[0])/baseW)+1):
        for i in range(len(baseP)):
            if primeB[i,j]:
                for i1 in range(len(baseP)):
                    j1=1-1*i1//(len(baseP)-1)+baseW*j*j+(baseP[i]+C1[i1,i])*j+np.dot(baseP[i],C1[i1,i])//baseW
                    primeB[i1,j1::(baseW*j+baseP[i])] =False
    return np.append([2,3,5],baseP[np.nonzero(primeB.T)[1][len(baseP):]]+baseW*np.nonzero(primeB.T)[0][len(baseP):])

N=1000000000
start = time.time()
primes = sieve(N)
    
end = time.time()

time_spent = end - start
print("\nPrimes count < ",N,": ", len(primes))
print("\nTime spent: ", time_spent, "seconds")

N=10^9 的结果：

Primes count <  1000000000 :  50847534

Time spent:  4.172985553741455 seconds

这是分段版本，占用内存更少：

import time
import numpy as np
from math import isqrt

def segmented_sieve(n , max_bW):
    dim_seg = 2**22
    bW = 30
    n_PB = 3
    n_PB_max = 5
    Primes_Base = np.array([2, 3, 5, 7, 11])
    if (max_bW > bW):
        while(n_PB < n_PB_max and (bW * Primes_Base[n_PB] <= max_bW) and n > bW**2):
            bW *= Primes_Base[n_PB]
            n_PB += 1
        Remainder_t = np.ones(bW, dtype=bool)
        for i in range(n_PB):
            Remainder_t[Primes_Base[i]::Primes_Base[i]] = False
        RW = np.empty(0, int)
        for i in range(2, bW):
            if (Remainder_t[i]):
                RW = np.append(RW, int(i))
        RW = np.append(RW, bW + 1)
        nR = len(RW)
    else:
        RW = np.array([7, 11, 13, 17, 19, 23, 29, 31])
        nR = 8
    k_end = N // bW + 1
    if (N % bW <= 1 and k_end > 0):
        k_end -= 1
    k_sqrt = isqrt(k_end // bW) + 2
    k_sqrt_sqrt = isqrt(k_sqrt // bW) + 2
    C_t = np.ones((nR, nR), dtype = int)
    RW_i = np.zeros((nR, nR), dtype = int)
    C_t = np.multiply(np.dot(np.diag(RW), C_t), np.dot(C_t, np.diag(RW)))
    C_t = C_t % bW
    C_t[C_t == 1] = C_t[C_t == 1] + bW
    for i in range(nR):
        RW_i[i,:] = RW[np.argwhere(C_t == RW[i])[:,1]]
    primeSqrt = np.ones((nR, k_sqrt), dtype = bool)
    for k in range(k_sqrt_sqrt):
        for i in range(nR):
            if primeSqrt[i, k]:
                for i1 in range(nR):
                    m = bW * k * k + (RW[i] + RW_i[i1,i]) * k + np.dot(RW[i], RW_i[i1,i]) // bW - i1 //(nR - 1)
                    primeSqrt[i1, m :: (bW * k + RW [i])] = False
    primes =  np.array(np.append(Primes_Base[:n_PB], RW[np.nonzero(primeSqrt.T)[1][::]] + bW * np.nonzero(primeSqrt.T)[0][::]))
    k_low = k_sqrt
    primeSeg = np.ones((nR, dim_seg), dtype = bool)
    while (k_low < k_end):
        for k in range(k_sqrt):
            for i in range(nR):
                p = bW * k + RW [i]
                if primeSqrt[i, k]:
                    for i1 in range(nR):
                        m = RW_i[i1,i] * k + np.dot(RW[i], RW_i[i1,i]) // bW - i1 // (nR - 1)
                        if (m < k_low):
                            m = p - ((k_low - m) % p);
                        else:
                            m -= k_low
                        if (m == p):
                            m = 0
                        if (m < dim_seg):
                            primeSeg[i1, m :: p] = False
        primes = np.append(primes, RW[np.nonzero(primeSeg.T)[1][::]] + bW * (np.nonzero(primeSeg.T)[0][::] + k_low))
        primeSeg.fill(True)        
        k_low += dim_seg
    return primes[primes < N]

N=1000000000
start = time.time()

primes = segmented_sieve(N, 30)
    
end = time.time()

time_spent = end - start
print("\nPrimes count < ",N,": ", len(primes))
print("\nTime spent: ", time_spent, "seconds")