记一次筛素数算法的优化
求出1..n中有多少个素数的算法我们通常用筛数法，一种比较简单的实现如下：
1 #include <stdio.h> 2 #include <stdlib.h> 3 4 int main(const int argc, const char *argv[], const char *envp[]) 5 { 6 int64_t n; 7 int64_t i, j; 8 int64_t *arr; 9 int64_t cnt = 0; 10 11 if (argc != 2) { 12 return -1; 13 } 14 15 sscanf(argv[1], "%lld", &n); 16 17 arr = (int64_t *)calloc(n + 1, sizeof(int64_t)); 18 19 for (i = 2; i <= n; i++) { 20 arr[i] = 1; 21 } 22 23 for (i = 2; i < n; i++) { // outer for-loop 24 25 if (arr[i] == 0) { 26 continue; 27 } 28 29 for (j = i + i; j <= n; j += i) { // inner for-loop 30 arr[j] = 0; 31 } 32 } 33 34 for (i = 2; i <= n; i++) { 35 cnt += arr[i] != 0; 36 } 37 38 printf("%lld/n", cnt); 39 40 return 0; 41 }
这里首先解释一下为什么外层循环不用取n，因为如果n是合数，显然n已经在i为n最小素因数的时候就已经被筛掉了；否则n是素数不用筛。
取n=100000000测测运行时间：
1 ~/ time ./test1 100000000 2 5761455 3 ./test1 100000000 2.92s user 0.23s system 99% cpu 3.157 total
现在我们注意到：每次进入内层for循环的时候，实际上对于素数i， 所有比i 2 小的i的倍数{j * i | j < i}都一定有一个比i更小的素因数 ，那我们其实可以将j从i 2 开始迭代，这样就可以节约一定的时间，简易实现如下：
1 #include <stdio.h> 2 #include <stdlib.h> 3 4 int main(const int argc, const char *argv[], const char *envp[]) 5 { 6 int64_t n; 7 int64_t i, j; 8 int64_t *arr; 9 int64_t cnt = 0; 10 11 if (argc != 2) { 12 return -1; 13 } 14 15 sscanf(argv[1], "%lld", &n); 16 17 arr = (int64_t *)calloc(n + 1, sizeof(int64_t)); 18 19 for (i = 2; i <= n; i++) { 20 arr[i] = 1; 21 } 22 23 for (i = 2; i < n; i++) { // outer for-loop 24 25 if (arr[i] == 0) { 26 continue; 27 } 28 29 for (j = i * i; j <= n; j += i) { // inner for-loop 30 arr[j] = 0; 31 } 32 } 33 34 for (i = 2; i <= n; i++) { 35 cnt += arr[i] != 0; 36 } 37 38 printf("%lld/n", cnt); 39 40 return 0; 41 }
仅仅改动了一个字符（'+' → '*'），有多大的时间收益呢？
1 ~/ time ./test2 100000000 2 5761455 3 ./test2 100000000 2.02s user 0.25s system 99% cpu 2.263 total
时间减少了近1s，几乎节约了1/3的时间，这个收效很明显啊。而且从输出来看，确实应该是正确的。（我知道单个测试样例不能证明程序的正确性，但是对于筛素数的算法而言，n取得足够大还能对，这就是一种佐证）
那我们再来思考思考还能不能再有所提升？
我们继续考虑内层循环，其对arr[j]的访问实际上相当没有空间连续性，并且因此也失掉了向量化的可能性。考虑到进入内层循环时所有比i小的素数的倍数都已经被筛过了，所以对于i的在[i, n/i]范围内的倍数中的也有很大部分被筛过了，也就是说 我们仅需要筛掉[i, n/i]中尚未被筛掉的数的i倍的那个数{j * i | j ∈ [i, n/i]，且j未被筛掉} 。这样的话我们首先可以得出一个推论：外层循环只需要取到sqrt(n)即可，实际中我直接用i * i <= n来作为外层循环的条件。另外，我们应当尽量消除条件分支，由于判断j是否被筛掉会造成一个if语句，所以我们可以改成内存循环中始终将arr[k]置0，若j未被筛掉则k = j * i，否则k = 0，综合起来就是k = arr[j] * j * i，这样就去掉了if语句并且利用了arr[0]留空这点，同时还使得如果不需要筛掉j * i的情况有很好的空间局部性。更进一步的，由于总是会计算arr[j] * j，所以我们的arr[]不要仅用来存放{0, 1}，而是直接用arr[j]存放j，这样k就是arr[j] * i。这样程序就变成了：
1 #include <stdio.h> 2 #include <stdlib.h> 3 4 int main(const int argc, const char *argv[], const char *envp[]) 5 { 6 int64_t n; 7 int64_t i, j; 8 int64_t *arr; 9 int64_t cnt = 0L; 10 11 if (argc != 2) { 12 return -1; 13 } 14 15 sscanf(argv[1], "%lld", &n); 16 17 arr = (int64_t *)calloc(n + 1, sizeof(int64_t)); 18 19 for (i = 2; i <= n; i++) { 20 arr[i] = i; 21 } 22 23 for (i = 2; i * i <= n; i++) { // outer for-loop 24 25 if (arr[i] == 0) { 26 continue; 27 } 28 29 for (j = n / i; j >= i; j--) { // inner for-loop 30 const int64_t k = arr[j] * i; 31 arr[k] = 0; 32 } 33 } 34 35 for (i = 2; i <= n; i++) { 36 cnt += arr[i] != 0; 37 } 38 39 printf("%lld/n", cnt); 40 41 return 0; 42 }
之所以这里内存循环变成了降序的迭代，是因为避免j和j * i同时都没被筛掉过的话，且j * i 2 < n的话，会引起错误。
这次程序改动较大，那有没有较大的时间上的收益呢？
1 ~/ time ./test3 100000000 2 5761455 3 ./test3 100000000 1.01s user 0.25s system 99% cpu 1.263 total
又减少了1s，只不过这1s优化来得比较麻烦。
正文到此结束
热门推荐
相关文章
近期评论
-
666
-
文章和留言都翻到11页了 没有OOM
-
-
朋友，翻页到11页，及以后，会出现OOM，无法访问
-
-
-
-
版本号是多少，你可以下载哪个代码仓库，jdk选1.8 直接跑就行
-
-
Loading...
![[HBLOG]公众号](https://www.liuhaihua.cn/img/qrcode_gzh.jpg)
