Eratosthenes的筛子 – find首席Python
只是澄清,这不是一个功课问题:)
我想为我正在构build的math应用程序find素数,并遇到了Eratosthenes的Sieve方法。
我已经用Python写了一个实现。 但速度非常慢。 比方说,如果我想find所有不到200万的素数。 它需要> 20分钟。 (我在此停止)。 我怎样才能加快速度呢?
def primes_sieve(limit): limitn = limit+1 primes = range(2, limitn) for i in primes: factors = range(i, limitn, i) for f in factors[1:]: if f in primes: primes.remove(f) return primes print primes_sieve(2000)
更新:我结束了对这个代码进行分析,发现花费了相当多的时间从列表中删除一个元素。 相当可以理解,考虑到它必须遍历整个列表(最坏的情况)来find元素,然后删除它,然后重新调整列表(也许一些副本呢?)。 无论如何,我抽出了字典的名单。 我的新实现 –
def primes_sieve1(limit): limitn = limit+1 primes = dict() for i in range(2, limitn): primes[i] = True for i in primes: factors = range(i,limitn, i) for f in factors[1:]: primes[f] = False return [i for i in primes if primes[i]==True] print primes_sieve1(2000000)
你不是很执行正确的algorithm:
在你的第一个例子中, primes_sieve
并不维护罢工/取消设置(如algorithm中)的素性标志列表,而是连续调整整数列表,这非常昂贵:从列表中删除项目需要将所有后续项目下降一个。
在第二个例子中, primes_sieve1
维护一个有素性标志的字典 ,这是一个正确方向的步骤,但它以不确定的顺序在字典上迭代,冗余地去除了因素中的因素(而不仅仅是素数的因素)algorithm)。 你可以通过对键进行sorting来解决这个问题,并且跳过非质数(已经使其速度提高了一个数量级),但是直接使用列表还是更有效率的。
正确的algorithm(用列表而不是字典)看起来像这样:
def primes_sieve2(limit): a = [True] * limit # Initialize the primality list a[0] = a[1] = False for (i, isprime) in enumerate(a): if isprime: yield i for n in xrange(i*i, limit, i): # Mark factors non-prime a[n] = False
(请注意,这也包括在素数的平方( i*i
)处开始非素数标记的algorithm优化,而不是它的double)。
def eratosthenes(n): multiples = [] for i in range(2, n+1): if i not in multiples: print (i) for j in range(i*i, n+1, i): multiples.append(j) eratosthenes(100)
从数组的开头(列表)移除需要移动所有的项目。 这意味着以这种方式从前面开始从列表中删除每个元素是一个O(n ^ 2)操作。
你可以更有效地做到这一点:
def primes_sieve(limit): limitn = limit+1 not_prime = set() primes = [] for i in range(2, limitn): if i in not_prime: continue for f in range(i*2, limitn, i): not_prime.add(f) primes.append(i) return primes print primes_sieve(1000000)
…或者,避免重新排列列表:
def primes_sieve(limit): limitn = limit+1 not_prime = [False] * limitn primes = [] for i in range(2, limitn): if not_prime[i]: continue for f in xrange(i*2, limitn, i): not_prime[f] = True primes.append(i) return primes
我意识到这并不是真正回答如何快速生成素数的问题,但也许有些人会发现这个select有趣:因为python提供了通过生成器的懒惰评估,eratosthenes的筛选可以完全按照说明实现:
def intsfrom(n): while True: yield n n += 1 def sieve(ilist): p = next(ilist) yield p for q in sieve(n for n in ilist if n%p != 0): yield q try: for p in sieve(intsfrom(2)): print p, print '' except RuntimeError as e: print e
try块在那里,因为algorithm运行,直到它吹到堆栈,没有try块显示回溯显示推实际输出你想看到屏幕。
通过结合许多爱好者的贡献(包括Glenn Maynard和MrHIDEn的评论),我在python 2中提出了以下代码片段:
def simpleSieve(sieveSize): #creating Sieve. sieve = [True] * (sieveSize+1) # 0 and 1 are not considered prime. sieve[0] = False sieve[1] = False for i in xrange(2,int(math.sqrt(sieveSize))+1): if sieve[i] == False: continue for pointer in xrange(i**2, sieveSize+1, i): sieve[pointer] = False # Sieve is left with prime numbers == True primes = [] for i in xrange(sieveSize+1): if sieve[i] == True: primes.append(i) return primes sieveSize = input() primes = simpleSieve(sieveSize)
在我的机器上为不同input的10次方input所需的时间是:
- 3:0.3ms
- 4:2.4毫秒
- 5:23毫秒
- 6:0.26秒
- 7:3.1 s
- 8:33秒
一个简单的速度黑客:当你定义variables“素数”,设置步骤2自动跳过所有的偶数,并设置起点为1。
那么你可以进一步优化,而不是我在素数中,用于我素数[:round(len(primes)** 0.5)]。 这将大大提高性能。 另外,可以删除以5结尾的数字,以进一步提高速度。
快多了:
def get_primes(n): m = n+1 numbers = [True for i in range(m)] for i in range(2, int(math.sqrt(n))): if numbers[i]: for j in range(i*i, m, i): numbers[j] = False primes = [] for i in range(2, m): if numbers[i]: primes.append(i) return primes start = time.time() primes = get_primes(10000) print(time.time() - start) print(get_primes(100))
使用filter
方法筛选数字列表,可以完成以下操作。
from math import sqrt def eratosthenes(limit): lst = range(1, limit) for i in range(2, int(sqrt(limit)) + 1): lst = filter(lambda x: x == i or x % i, lst) # sieve return lst print eratosthenes(2000000)
我的实现:
import math n = 100 marked = {} for i in range(2, int(math.sqrt(n))): if not marked.get(i): for x in range(i * i, n, i): marked[x] = True for i in range(2, n): if not marked.get(i): print i
这是一个更有记忆效率的版本(和:适当的筛选,而不是审判部门)。 基本上,不是保留所有数字的数组,而是去掉那些不是最主要的数字,这样就保留了一系列的计数器 – 每个发现的数字都有一个 – 并且在假定的素数之前跳过它们。 这样,它使用与素数成比例的存储,而不是达到最高素数。
import itertools def primes(): class counter: def __init__ (this, n): this.n, this.current, this.isVirgin = n, n*n, True # isVirgin means it's never been incremented def advancePast (this, n): # return true if the counter advanced if this.current > n: if this.isVirgin: raise StopIteration # if this is virgin, then so will be all the subsequent counters. Don't need to iterate further. return False this.current += this.n # pre: this.current == n; post: this.current > n. this.isVirgin = False # when it's gone, it's gone return True yield 1 multiples = [] for n in itertools.count(2): isPrime = True for p in (m.advancePast(n) for m in multiples): if p: isPrime = False if isPrime: yield n multiples.append (counter (n))
你会注意到primes()
是一个生成器,所以你可以把结果保存在一个列表中,或者你可以直接使用它们。 这是前n
素数:
import itertools for k in itertools.islice (primes(), n): print (k)
而且,为了完整性,这里有一个计时器来衡量性能:
import time def timer (): t, k = time.process_time(), 10 for p in primes(): if p>k: print (time.process_time()-t, " to ", p, "\n") k *= 10 if k>100000: return
万一你想知道,我也写了primes()
作为一个简单的迭代器(使用__iter__
和__next__
),它以几乎相同的速度运行。 我也感到惊讶!
因为速度我更喜欢NumPy。
import numpy as np # Find all prime numbers using Sieve of Eratosthenes def get_primes1(n): m = int(np.sqrt(n)) is_prime = np.ones(n, dtype=bool) is_prime[:2] = False # 0 and 1 are not primes for i in range(2, m): if is_prime[i] == False: continue is_prime[i*i::i] = False return np.nonzero(is_prime)[0] # Find all prime numbers using brute-force. def isprime(n): ''' Check if integer n is a prime ''' n = abs(int(n)) # n is a positive integer if n < 2: # 0 and 1 are not primes return False if n == 2: # 2 is the only even prime number return True if not n & 1: # all other even numbers are not primes return False # Range starts with 3 and only needs to go up the square root # of n for all odd numbers for x in range(3, int(n**0.5)+1, 2): if n % x == 0: return False return True # To apply a function to a numpy array, one have to vectorize the function def get_primes2(n): vectorized_isprime = np.vectorize(isprime) a = np.arange(n) return a[vectorized_isprime(a)]
检查输出:
n = 100 print(get_primes1(n)) print(get_primes2(n)) [ 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] [ 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]
比较Eratosthenes Sieve的速度和Jupyter Notebook上的蛮力。 Eratosthenes筛比539倍的速度比百万元的蛮力。
%timeit get_primes1(1000000) %timeit get_primes2(1000000) 4.79 ms ± 90.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each) 2.58 s ± 31.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Here is my c++ solution #include <iostream> #include <cmath> #include <iomanip> using namespace std; const int MAX_NUMBER = 1000;// Max size of sequence of numbers. const int NUMBERS_PER_LINE = 20; // Max sequence of numbers printed per line in output. /** Initializes input boolean array to default value true. @param input_arr The array to be initialized. */ void arr_initialzer(bool[]); /** Sets Composite Number indexes to false, while Prime Number indexes are left true. Calculations are done by using Sieve of Eratosthenes Algorithm. @param input_arr The Input array for computing prime numbers. */ void compute_primes(bool []); /** Prints the prime number indexes. @param input_arr The Input array passed for printing. */ void print_primes(bool[]); int main() { bool seq_arr[MAX_NUMBER +1]; //Initialize the seq_arr to default value true. arr_initialzer(seq_arr); // Compute prime numbers compute_primes(seq_arr); // Print prime numbers print_primes(seq_arr); return 0; } void arr_initialzer(bool input_arr[]) { for (int i = 2;i < MAX_NUMBER+1;i++) { input_arr[i] = true; } } void compute_primes(bool input_arr[]) { int k = 1; // Multiplying factor. for (int i = 2; i < sqrt(MAX_NUMBER); i++) { if (input_arr[i]) { for (int j = i*i; j < MAX_NUMBER +1; j = i*i + i*k++) { input_arr[j] = false; } k = 1; } } } void print_primes(bool input_arr[]) { int check_number_count = 0; for (int i = 2;i < MAX_NUMBER+1;i++) { if (input_arr[i]) { if (check_number_count == MAX_NUMBER) { cout << endl; check_number_count = 0; } cout <<setw(4)<<i; check_number_count++; } } }