素数筛法系列之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取代.