素数筛法系列之1 简单筛法实现

    让你写一个程序, 计算10亿(10^9)内素数(写入内存, 不必输出), 你能在2G的主流cpu上10秒做到?

.看看这个作者的实现, 其实0.5秒就够了. http://code.google.com/p/primesieve/
 

    如何快速枚举n以内的素数, 目前最快的算法还是 sieve of Eratosthenes

要高效的实现起来并不那么容易, 一定程度上可考察你的基本功, 尤其是如何去优化程序

方面的功底.

下面是我写的一端小程序, 看看有什么问题.

/*

    为节省内存使用char存相邻素数的差, 尤其是n比较大,并以后还有大量的遍历操作,

可以提高cache命中率, 但代码复杂度也高不少.

*/

static unsigend char Prime[10000];

/****

the classic sieve of Eratosthenes implementation by bit packing

all prime less than sqrt(n) will be saved in Prime[] by difference

Prime[1] = 2, Prime[2] = 3 - 2 = 1, Prime[3] = 5 - 3, ....

Prime[i] is the difference of the adjacent prime:

Prime[i] = (ith Prime) - ((i-1) th Prime);

*/

static int simpleEratoSieve1(const uint n)

{

    int primes = 1;

    uchar bitarray[(1000000 >> (3 + 1)) + 10] = {0}; // 注意stackoverflow!!!

    assert(n < sizeof(bitarray) * 16);

    uint lastprime = Prime[primes++] = 2;

    for (uint p = 3; p <= n; p += 2) {

        //bit position with vlaue 0 is prime number

        if (!TST_BIT(bitarray, p / 2)) {

            Prime[primes++] = p - lastprime;

            lastprime = p;

            if (p > maxp / p) {

                continue;

            }

            for (uint j = p * p / 2; j <= maxp / 2; j += p) {

                SET_BIT(bitarray, j);

            }

        }

    }

    //pack the last two number

    Prime[primes + 0] = Prime[primes + 1] = 255;

    return primes;

}

看完代码之后有几个问题.

1. 难道相邻素数差小于uchar最大值255?

2. 如果需要计算更大范围内素数怎么办?

3. 如果要直接存每一个素数而不是相邻差?

4. 为什么用位运算而不是bool数组实现?

上面的问题解答如下

1. 在10^9内只有三对相邻素数差值大于UCHAR_MAX = 255, 目前已知最大素数间隔不超过2000

解决办法是把Prime的定义修改为

    static unsigend short Prime[10000];

2. 修改内存bitarray实现

        uchar* bitarray = new uchar [(n >> (3 + 1)) + 9];

        memset(bitarray, 0, (n >> (3 + 1)) + 9);

3. 这个更简单了,去除lastprime

4. 较大范围内存性能更好(cache friendly), 不信你可以试试计算n = 1e8或更大

    

static int simpleEratoSieve2(const uint64 n)

{

    int primes = 1;

    uchar* bitarray = new uchar [(n >> (3 + 1)) + 9];

    memset(bitarray, 0, (n >> (3 + 1)) + 9);

    Prime[primes++] = 2;

    for (uint p = 3; p <= n; p += 2) {

        //bit position with vlaue 0 is prime number

        if (!TST_BIT(bitarray, p / 2)) {

            Prime[primes++] = p;

            if (p > n / p) {

                continue;

            }

            for (uint j = p * p / 2; j <= n / 2; j += p) {

                SET_BIT(bitarray, j);

            }

        }

    }

    //pack the last two number

    Prime[primes + 0] = -1u;

    delete  bitarray;

    return primes;

}

simpleEratoSieve2简单的多,功能也强大, 能计算10^12以内素数吗? 自己去试试看, 会发现很多问题.

如果计算10^12以内素数, 只需先计算10^6以内素数表, simpleEratoSieve1构造小范围素数表

, 再用区间分段的算法去统计每一个区间的个数, 请看下回分解.

结论 simpleEratoSieve1已足够满足需求, simpleEratoSieve2会被the segmented sieve of Eratosthenes取代.