Sieve of Eratosthenes筛法

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秒左右

Sieve of Eratosthenes
Sieve of Eratosthenes

于是完整的算法如下,

//
//  main.cpp
//  sieve
//
//  Created by BlueCocoa on 16/1/27.
//  Copyright © 2016 BlueCocoa. 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;
}
声明: 本文为0xBBC原创, 转载注明出处喵~

《Sieve of Eratosthenes筛法》有1个想法

发表评论

电子邮件地址不会被公开。 必填项已用*标注