<p>我学到了一个技巧<a href="https://stackoverflow.com/q/62897329/5987">just yesterday</a>——如果你把数字分成6组,6组中只有2个可能是素数。其他的可以被2或3平分。这意味着跟踪6个数字的素性只需要2位;一个包含8位的字节可以跟踪24个数字的素数!这大大减少了筛的内存需求</p>
<p>在Windows10上的Python3.7.5 64位中,以下代码没有超过36.4MB</p>
<pre><code>remainder_bit = [0, 0x01, 0, 0, 0, 0x02,
0, 0x04, 0, 0, 0, 0x08,
0, 0x10, 0, 0, 0, 0x20,
0, 0x40, 0, 0, 0, 0x80]
def is_prime(xs, a):
if a <= 3:
return a > 1
index, rem = divmod(a, 24)
bit = remainder_bit[rem]
if not bit:
return False
return not (xs[index] & bit)
def sieve_of_eratosthenes(xs, n):
count = (n // 3) + 1 # subtract out 1 and 4, add 2 3 and 5
p = 5
while p*p <= n:
if is_prime(xs, p):
for i in range(5 * p, n + 1, p):
index, rem = divmod(i, 24)
bit = remainder_bit[rem]
if bit and not (xs[index] & bit):
xs[index] |= bit
count -= 1
p += 2
if is_prime(xs, p):
for i in range(5 * p, n + 1, p):
index, rem = divmod(i, 24)
bit = remainder_bit[rem]
if bit and not (xs[index] & bit):
xs[index] |= bit
count -= 1
p += 4
return count
def init_sieve(n):
return bytearray((n + 23) // 24)
n = 100000000
xs = init_sieve(n)
sieve_of_eratosthenes(xs, n)
5761455
sum(is_prime(xs, i) for i in range(n+1))
5761455
</code></pre>
<hr/>
编辑:了解其工作原理的关键是,筛子会创建重复模式。对于素数2和3,模式每2*3或6个数字重复一次,在这6个数字中,4个已经被渲染为不可能是素数,只剩下2个。没有什么限制你选择素数来生成模式,也许除了收益递减定律。我决定尝试在混音中加入5,使图案每2*3*5=30个数字重复一次。在这30个数字中,只有8个可以是素数,这意味着每个字节可以跟踪30个数字,而不是上面的24个!这使您在内存使用方面有20%的优势。
<p>这是更新后的代码。我还简化了一点,去掉了素数的计数</p>
<pre><code>remainder_bit30 = [0, 0x01, 0, 0, 0, 0, 0, 0x02, 0, 0,
0, 0x04, 0, 0x08, 0, 0, 0, 0x10, 0, 0x20,
0, 0, 0, 0x40, 0, 0, 0, 0, 0, 0x80]
def is_prime(xs, a):
if a <= 5:
return (a > 1) and (a != 4)
index, rem = divmod(a, 30)
bit = remainder_bit30[rem]
return (bit != 0) and not (xs[index] & bit)
def sieve_of_eratosthenes(xs):
n = 30 * len(xs) - 1
p = 0
while p*p < n:
for offset in (1, 7, 11, 13, 17, 19, 23, 29):
p += offset
if is_prime(xs, p):
for i in range(p * p, n + 1, p):
index, rem = divmod(i, 30)
if index < len(xs):
bit = remainder_bit30[rem]
xs[index] |= bit
p -= offset
p += 30
def init_sieve(n):
b = bytearray((n + 30) // 30)
return b
</code></pre>