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;
}