用C++实现的增强Eratosthenes筛法程序

运行示例

PS H:\Read\num\x64\Release> .\eSievePro
Eratosthenes sieve: a method to find out all primes below the number that you specify here please: 100

Only the sum of all primes needed [y/n](y as default): n

Start to work out all primes below 100...
25 primes found in 0 milliseconds.

2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97
PS H:\Read\num\x64\Release> .\eSievePro
Eratosthenes sieve: a method to find out all primes below the number that you specify here please: 12345678

Only the sum of all primes needed [y/n](y as default):

Start to work out the sum of all primes below 12345678...
809227 primes found in 63 milliseconds.

PS H:\Read\num\x64\Release> .\eSievePro
Eratosthenes sieve: a method to find out all primes below the number that you specify here please: 123456789

Only the sum of all primes needed [y/n](y as default):

Start to work out the sum of all primes below 123456789...
7027260 primes found in 1031 milliseconds.

PS H:\Read\num\x64\Release> .\eSievePro
Eratosthenes sieve: a method to find out all primes below the number that you specify here please: 1234567890

Only the sum of all primes needed [y/n](y as default):

Start to work out the sum of all primes below 1234567890...
62106578 primes found in 12375 milliseconds.

PS H:\Read\num\x64\Release> .\eSievePro
Eratosthenes sieve: a method to find out all primes below the number that you specify here please: 1234567

Only the sum of all primes needed [y/n](y as default):

Start to work out the sum of all primes below 1234567...
95360 primes found in 0 milliseconds.

PS H:\Read\num\x64\Release> 

和用C++实现的Eratosthenes筛法程序对比，性能有显著提升。

 

增强Eratosthenes筛法的C++程序实现

为说明方便，把用C++实现的Eratosthenes筛法程序里的实现称为eSieve，而把这里的新实现称为eSievePro。

宏定义和全局量定义

typedef unsigned char u8;
typedef unsigned long ulong;
static ulong s_last = 0;
static u8* s_pOdd = NULL;
static std::vector<ulong> s_vecPrime;

对照eSieve的相应部分，这里只是把之前的s_pAll换成了s_pOdd。

main函数实现

int main()
{
    printf(" Eratosthenes sieve: a method to find out all primes below the number that you specify here please: ");
    std::string strInput;
    getline(std::cin, strInput);
    ulong raw_last = strtoul(strInput.c_str(), 0, 10);
    if (raw_last <= 2) {
        printf("\n Wrong input.\n");
        return 0;
    }
    printf("\n Only the sum of all primes needed [y/n](y as default): ");
    getline(std::cin, strInput);
    bool bDetail = (strInput == "n");
    if (bDetail)
        printf("\n Start to work out all primes below %u...\n", raw_last);
    else
        printf("\n Start to work out the sum of all primes below %u...\n", raw_last);
    if (!eSievePro(raw_last, bDetail))
        return 0;
    if (bDetail)
        showDetails();
    return 0;
}

这里的main函数实现保留了交互输入的实现代码，而把增强的Eratosthenes筛法实现放到了单独的eSievePro函数里。

eSievePro函数实现

bool eSievePro(ulong raw_last, bool bDetail)
{
    DWORD tickBegin = GetTickCount();
    s_last = raw_last / 2;
    s_pOdd = new u8[s_last];
    if (!s_pOdd) {
        printf("Lack of memory.\n");
        return false;
    }
    ulong sum = 1;
    ulong curPrime = 2;
    if (bDetail)
        s_vecPrime.push_back(curPrime);
    if (raw_last == 3)
        goto Ending;

    memset(s_pOdd, 1, s_last);
    ulong curPrimeIdx = 0;
    while (true) {
        renewCurrentPrime(curPrimeIdx);
        ++sum;
        curPrime = (curPrimeIdx + curPrimeIdx) + 1;
        if (bDetail)
            s_vecPrime.push_back(curPrime);
        s_pOdd[curPrimeIdx - 1] = 0;
        ulong sqr = curPrime * curPrime;
        if (sqr > raw_last)
            break;
        ulong step = curPrime + curPrime;
        for (ulong idx = sqr; idx < raw_last; idx += step) {
            s_pOdd[(idx - 1) / 2 - 1] = 0;
        }
    }
    /// pick up all the left primes
    for (ulong idx = curPrime + 2; idx < raw_last; idx += 2) {
        ulong halfIdx = (idx - 1) / 2 - 1;
        if (s_pOdd[halfIdx] == 1) {
            ++sum;
            if (bDetail)
                s_vecPrime.push_back(idx);
        }
    }

Ending:
    printf(" %u primes found in %u milliseconds.\n\n", sum, GetTickCount() - tickBegin);
    delete []s_pOdd;
    return true;
}

由

 4     s_last = raw_last / 2;
 5     s_pOdd = new u8[s_last];

这两行代码可以看出，s_pOdd的动态申请的空间是原来的一半，就是说原来的动态数组s_pAll记录的是1、2、3、……、upperLimit（upperLimit指代交互输入的上界值）所有自然数的素合状态（1表示素数，0表示合数）；而这里的动态数组s_pOdd只用一半的空间来记录1到upperLimit之内的所有奇数的素合状态。这样处理不但节省了存储空间，而且直接省去了Eratosthenes筛法中用2筛去所有偶数的处理过程。

如下这几行筛法实现代码

26         ulong sqr = curPrime * curPrime;

29         ulong step = curPrime + curPrime;
30         for (ulong idx = sqr; idx < raw_last; idx += step) {
31             s_pOdd[(idx - 1) / 2 - 1] = 0;
32         }

相比之前的实现代码

36         for (int idx = curPrime - 1; idx <= s_last - 1; idx += curPrime) {
37             s_pAll[idx] = 0;
38         }

有两个改进之处。用一个实例说明一下：比如用Eratosthenes筛法筛去1到1000之内的合数，依次筛去2的倍数、3的倍数、5的倍数、7的倍数，当开始筛去11的倍数时，会依次筛去22、33、44、55、66、77、88、99、110、121、132、143、……，但是22、33、44、55、66、77、88、99、110这些数因为小于11*11，必然有小于11的素因素，说明这些数必然会在此前筛去2、3、5、7等素数的倍数的过程中被筛掉，因此当筛去一个素数的倍数时，从该素数的平方为起点开始筛会减少重复，在上面的代码对比上用标黄部分表示。

另一个改进之处，在上面的代码对比上用标蓝部分表示。仍沿用刚才的实例说明，当从121开始筛去11的倍数时，会依次筛去121、132、143、154、165、176、187、……，从这个列表可以看出奇数偶数交替出现，这也很好解释，因为筛法当前用到的素数（curPrime）是一个奇数，设为2n+1，它的平方也是奇数（即4nn+4n+1），再加上2n+1就会变成偶数（即4nn+6n+2），再加上2n+1又会变成奇数（4nn+8n+3），……。因为这些偶数在最早筛去2的倍数时就已经被筛去了，而且2之后的素数都是奇数，所以可以把2之后的筛法步长调整为当时使用素数的两倍，这样性能可以提升一倍。

当然，上述两个改进之处并不能彻底避免多次重复筛去一个合数的问题。比如，刚才的实例中使用步长22会依次筛去121、143、165、187，但是里面的165=11*5*3，在之前已经被3筛去过一次，又被5筛去过一次。

renewCurrentPrime和showDetails

bool renewCurrentPrime(ulong& primeIdx)
{
    while (primeIdx < s_last) {
        ++primeIdx;
        if (s_pOdd[primeIdx - 1] == 1)
            return true;
    }
    return false;
}

这里的renewCurrentPrime和eSieve实现的renewCurrentPrime几乎一样，只是因为s_last和输入参数的含义有变化以及s_pOdd与此前的的s_pAll的不同含义而有些微不同。

showDetails的实现则完全没有变动，这里不再重复放了。

其他

https://github.com/readalps/eSievePro上放了eSievePro实现的源码文件，以及一个运行结果文件。