没有 2 和 3 的倍数的埃拉托色尼筛法
Sieve of Eratosthenes without multiples of 2 and 3
我正在尝试实现本文 Faster Sieve of Eratosthenes 中描述的算法。我很理解这个想法,但我无法理解如何通过 python 代码实现它。
经过一些工作后,我找到了一种将筛子中的索引转换为数字本身的方法:
number = lambda i: 3 * (i + 2) - 1 - (i + 2) % 2
但主要问题是我在获得质数后必须做的跳跃。文章解释为:
6np ± p, where p is prime and n - some natural number.
有没有一种方法可以使用最后找到的素数的索引来描述这种想法的跳跃?
提前致谢。
P.S。
有 implementation in Objective-C
我是编程新手,只能看懂python和js代码
如果您了解 numpy 以及 Python,请查看 primesfrom2to
的实现,摘自 this answer in Whosebug。
def primesfrom2to(n):
#
""" Input n>=6, Returns a array of primes, 2 <= p < n """
sieve = np.ones(n/3 + (n%6==2), dtype=np.bool)
sieve[0] = False
for i in xrange(int(n**0.5)/3+1):
if sieve[i]:
k=3*i+1|1
sieve[ ((k*k)/3) ::2*k] = False
sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False
return np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]
在我链接到的那个答案中,这个例程是建立素数列表最快的。我在自己的探索素数的代码中使用了它的变体。详细解释这一点需要很多 space,但它构建了一个 sieve
,省略了 2 和 3 的倍数。只有两行代码(开始 sieve[
和结束 = False
) 标出一个新发现的素数的倍数。我认为这就是您所说的 "jumps... after getting prime" 的意思。代码很棘手,但正在努力学习。该代码适用于 Python 2,遗留 Python.
这里是我自己的一些Python3代码,注释比较多。大家可以拿来对比一下。
def primesieve3rd(n):
"""Return 'a third of a prime sieve' for values less than n that
are in the form 6k-1 or 6k+1.
RETURN: A numpy boolean array. A True value means the associated
number is prime; False means 1 or composite.
NOTES: 1. If j is in that form then the index for its primality
test is j // 3.
2. The numbers 2 and 3 are not handled in the sieve.
3. This code is based on primesfrom2to in
<
fastest-way-to-list-all-primes-below-n-in-python/
3035188#3035188>
"""
if n <= 1: # values returning an empty array
return np.ones(0, dtype=np.bool_)
sieve = np.ones(n // 3 + (n % 6 == 2), dtype=np.bool_) # all True
sieve[0] = False # note 1 is not prime
for i in range(1, (math.ceil(math.sqrt(n)) + 1) // 3): # sometimes large
if sieve[i]:
k = 3 * i + 1 | 1 # the associated number for this cell
# Mark some of the stored multiples of the number as composite
sieve[k * k // 3 :: 2 * k] = False
# Mark the remaining stored multiples (k times next possible prime)
sieve[k * (k + 4 - 2*(i&1)) // 3 :: 2 * k] = False
return sieve
def primesfrom2to(n, sieve=None):
"""Return an array of prime numbers less than n.
RETURN: A numpty int64 (indices type) array.
NOTES: 1. This code is based on primesfrom2to in
<
fastest-way-to-list-all-primes-below-n-in-python/
3035188#3035188>
"""
if n <= 5:
return np.array([2, 3], dtype=np.intp)[:max(n - 2, 0)]
if sieve is None:
sieve = primesieve3rd(n)
elif n >= 3 * len(sieve) + 1 | 1: # the next number to note in the sieve
raise ValueError('Number out of range of the sieve in '
'primesfrom2to')
return np.r_[2, 3, 3 * np.nonzero(sieve)[0] + 1 | 1]
有什么不明白的就问我吧
我正在尝试实现本文 Faster Sieve of Eratosthenes 中描述的算法。我很理解这个想法,但我无法理解如何通过 python 代码实现它。
经过一些工作后,我找到了一种将筛子中的索引转换为数字本身的方法:
number = lambda i: 3 * (i + 2) - 1 - (i + 2) % 2
但主要问题是我在获得质数后必须做的跳跃。文章解释为:
6np ± p, where p is prime and n - some natural number.
有没有一种方法可以使用最后找到的素数的索引来描述这种想法的跳跃?
提前致谢。
P.S。 有 implementation in Objective-C 我是编程新手,只能看懂python和js代码
如果您了解 numpy 以及 Python,请查看 primesfrom2to
的实现,摘自 this answer in Whosebug。
def primesfrom2to(n):
#
""" Input n>=6, Returns a array of primes, 2 <= p < n """
sieve = np.ones(n/3 + (n%6==2), dtype=np.bool)
sieve[0] = False
for i in xrange(int(n**0.5)/3+1):
if sieve[i]:
k=3*i+1|1
sieve[ ((k*k)/3) ::2*k] = False
sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False
return np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]
在我链接到的那个答案中,这个例程是建立素数列表最快的。我在自己的探索素数的代码中使用了它的变体。详细解释这一点需要很多 space,但它构建了一个 sieve
,省略了 2 和 3 的倍数。只有两行代码(开始 sieve[
和结束 = False
) 标出一个新发现的素数的倍数。我认为这就是您所说的 "jumps... after getting prime" 的意思。代码很棘手,但正在努力学习。该代码适用于 Python 2,遗留 Python.
这里是我自己的一些Python3代码,注释比较多。大家可以拿来对比一下。
def primesieve3rd(n):
"""Return 'a third of a prime sieve' for values less than n that
are in the form 6k-1 or 6k+1.
RETURN: A numpy boolean array. A True value means the associated
number is prime; False means 1 or composite.
NOTES: 1. If j is in that form then the index for its primality
test is j // 3.
2. The numbers 2 and 3 are not handled in the sieve.
3. This code is based on primesfrom2to in
<
fastest-way-to-list-all-primes-below-n-in-python/
3035188#3035188>
"""
if n <= 1: # values returning an empty array
return np.ones(0, dtype=np.bool_)
sieve = np.ones(n // 3 + (n % 6 == 2), dtype=np.bool_) # all True
sieve[0] = False # note 1 is not prime
for i in range(1, (math.ceil(math.sqrt(n)) + 1) // 3): # sometimes large
if sieve[i]:
k = 3 * i + 1 | 1 # the associated number for this cell
# Mark some of the stored multiples of the number as composite
sieve[k * k // 3 :: 2 * k] = False
# Mark the remaining stored multiples (k times next possible prime)
sieve[k * (k + 4 - 2*(i&1)) // 3 :: 2 * k] = False
return sieve
def primesfrom2to(n, sieve=None):
"""Return an array of prime numbers less than n.
RETURN: A numpty int64 (indices type) array.
NOTES: 1. This code is based on primesfrom2to in
<
fastest-way-to-list-all-primes-below-n-in-python/
3035188#3035188>
"""
if n <= 5:
return np.array([2, 3], dtype=np.intp)[:max(n - 2, 0)]
if sieve is None:
sieve = primesieve3rd(n)
elif n >= 3 * len(sieve) + 1 | 1: # the next number to note in the sieve
raise ValueError('Number out of range of the sieve in '
'primesfrom2to')
return np.r_[2, 3, 3 * np.nonzero(sieve)[0] + 1 | 1]
有什么不明白的就问我吧