Sieve of Eratosthenes筛法是一种找质数的算法。我们知道,对于一个合数,那么它必定可以分解为两个质数相乘,比如21=3x7。那么我们也可以理解为,一个合数必定是一个质数的倍数,对于刚才的例子,21就是质数3的7倍(当然也可以说是质数7的3倍)。
那么我们求给定范围内的质数就可以利用这一点——先把所有的数列出来,然后从2开始,把所有数的倍数都去掉,这样我们就能留下给定范围内的质数了。举个例子,求42以内的所有质数,那么我们的算法大致上是这样操作的:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ... 42 1 2 3 0 5 0 7 0 9 0 11 0 13 0 15 0 ... 42 1 2 3 0 5 0 7 0 0 0 11 0 13 0 0 0 ... 42 ...
我们在第一行列出了给定范围内的所有数字,1到42
在第二行里,我们去掉了2的所有倍数,在算法中我们是将它们标记为0
接下来要做的事情和第二步几乎一样,每次都找出一个最小的没有被标记为0的数,去掉它的倍数(这个数小于等于给定的范围)即可。
实际上,我们不必从2循环到42,只需要循环到floor(sqrt(42))即可。因为floor(sqrt(42)) = 6,且6x7=42,如果继续开始循环的话,我们将得到7x2=14, 7x3=21, 7x4=28, 7x5=35, 7x6=42,可以看出这些数字都是我们在循环2, 3, 5的时候,已经被标记为0的。
并且,从上面的分析中我们可以看出,对于每次循环到的数(以下称为i),我们从i*i开始标记即可。因为对于任一满足1<=t<=i的数 i*(i-t) 来说,都已经在我们先前循环到(i-t)时被标记过了。举个例子,对于5而言,我们应该从5x5=25开始标记,因为5x2=10(在i=2时被标记),5x3=15(在i=3时被标记),5x4=20(在i=2时被标记)。
再优化一步,因为我们只需要标记一个数是否是另一个数的倍数,所以我们不需要用int []来保存,也不需要用bool []来保存(bool类型占用的空间和char一样,1个字节),我们将直接用bit来标记。
首先是算出给定范围n所需要的位长
// 计算需要多少个字节, (n + 8)确保有足够的位 size_t length = (n + 8) >> CHAR_BIT_LOG;
接下来就是bit标记,i从2到floor(sqrt(n)),只要i对应的位没有被标记,就从i*i开始标记。
#define CHAR_BIT_LOG 3 #define MASK (~(~0 << CHAR_BIT_LOG)) #define setbit(a, x) ((a)[(x) >> CHAR_BIT_LOG] |= 1 << ((x) & MASK)) #define isset(a, x) ((a)[(x) >> CHAR_BIT_LOG] & (1 << ((x) & MASK))) for (int i = 2; i <= (int)floor(sqrt(n)); ++i) { if (!isset(result, i)) { int j = i*i; while (j <= n) { setbit(result, j); j += i; } } }
当n为400000000时,Sieve筛选法耗时3.56秒左右
于是完整的算法如下,
// // main.cpp // sieve // // Created by Ryza 16/1/27. // Copyright © 2016[data deleted]. All rights reserved. // #include <limits.h> #include <iostream> #include <functional> #include <chrono> #include <cmath> using namespace std; #define CHAR_BIT_LOG 3 #define MASK (~(~0 << CHAR_BIT_LOG)) #define setbit(a, x) ((a)[(x) >> CHAR_BIT_LOG] |= 1 << ((x) & MASK)) #define isset(a, x) ((a)[(x) >> CHAR_BIT_LOG] & (1 << ((x) & MASK))) char * sieve(unsigned long long n) { char * result = NULL; if (n >= 2) { // 计算需要多少个字节, (n + 8)确保有足够的位 size_t length = (n + 8) >> 3; result = (char *)malloc(length); memset(result, 0, length); for (int i = 2; i <= (int)floor(sqrt(n)); ++i) { if (!isset(result, i)) { int j = i*i; while (j <= n) { setbit(result, j); j += i; } } } } return result; } template <class Return, class Param> auto runtime(function<Return(Param)> func, Param param) -> Return { Return result; auto start = chrono::high_resolution_clock::now(); result = func(param); auto end = chrono::high_resolution_clock::now(); printf("%lf\n", ((chrono::duration<double>)((end - start))).count()); return result; } int main(int argc, const char * argv[]) { int upperbound = 42; char * result = runtime<char *, unsigned long long>(sieve, upperbound); for (int i = 2; i <= upperbound; i++) if (!isset(result, i)) printf("%d ", i); printf("\n"); free(result); return 0; }