C++阿特金筛法忽略了一些素数

Question

最近，我一直在研究一个使用阿特金筛法的 C++ 素数生成器（http://en. wikipedia.org/wiki/Sieve_of_atkin ）来生成它的素数。我的目标是能够生成任何 32 位数字。我将主要用它来解决项目欧拉问题。大多数情况下这只是一个夏季项目。

该程序使用位板来存储素数：即一系列 1 和 0，例如第 11 位为 1，第 12 位为 0，第 13 位为 1 等。为了有效使用内存，这实际上是和字符数组，每个字符包含 8 位。我使用标志和按位运算符来设置和检索位。该算法的原理很简单：使用一些我不假装理解的方程进行第一次传递，以定义一个数字是否被视为“素数”。这在很大程度上会得到正确的答案，但一些非质数将被标记为质数。因此，当迭代列表时，您将刚刚找到的素数的所有倍数设置为“非素数”。这样做有一个方便的优点，即素数越大，所需的处理器时间就越少。

我已经完成了 90%，但有一个问题： 一些素数丢失了。

通过检查位板，我确定这些素数在第一次传递期间被省略，这基本上为许多方程的每个解决方案切换了一个数字（请参阅维基百科条目）。我一次又一次地检查了这段代码。我什至尝试增加维基百科文章中显示的范围，这效率较低，但我认为可能会触及一些我以某种方式省略的数字。没有任何效果。这些数字只是评估为不是质数。我的大部分测试都是针对 128 以下的所有素数。在这个范围内，这些是被省略的素数：

23 和 59。

我毫不怀疑，在更高的范围内，会丢失更多素数（只是不想计算）通过他们所有人）。我不知道为什么这些都消失了，但它们确实存在。这两个素数有什么特别之处吗？我已经进行了两次和三次检查，发现并修复了错误，但它仍然可能是我错过的一些愚蠢的东西。

无论如何，这是我的代码：

#include <iostream>
#include <limits.h>
#include <math.h>

using namespace std;

const unsigned short DWORD_BITS = 8;

unsigned char flag(const unsigned char);
void printBinary(unsigned char);


class PrimeGen
{
    public:
        unsigned char* sieve;
        unsigned sievelen;
        unsigned limit;
        unsigned bookmark;


        PrimeGen(const unsigned);

        void firstPass();
        unsigned next();

        bool getBit(const unsigned);
        void onBit(const unsigned);
        void offBit(const unsigned);
        void switchBit(const unsigned);

        void printBoard();
};


PrimeGen::PrimeGen(const unsigned max_num)
{
    limit = max_num;
    sievelen = limit / DWORD_BITS + 1;
    bookmark = 0;

    sieve = (unsigned char*) malloc(sievelen);
    for (unsigned i = 0; i < sievelen; i++) {sieve[i] = 0;}

    firstPass();
}


inline bool PrimeGen::getBit(const unsigned index)
{
    return sieve[index/DWORD_BITS] & flag(index%DWORD_BITS);
}


inline void PrimeGen::onBit(const unsigned index)
{
    sieve[index/DWORD_BITS] |= flag(index%DWORD_BITS);
}


inline void PrimeGen::offBit(const unsigned index)
{
    sieve[index/DWORD_BITS] &= ~flag(index%DWORD_BITS);
}


inline void PrimeGen::switchBit(const unsigned index)
{
    sieve[index/DWORD_BITS] ^= flag(index%DWORD_BITS);
}


void PrimeGen::firstPass()
{
    unsigned nmod,n,x,y,xroof, yroof;

    //n = 4x^2 + y^2
    xroof = (unsigned) sqrt(((double)(limit - 1)) / 4);
    for(x = 1; x <= xroof; x++){
        yroof = (unsigned) sqrt((double)(limit - 4 * x * x));
        for(y = 1; y <= yroof; y++){
            n = (4 * x * x) + (y * y);
            nmod = n % 12;
            if (nmod == 1 || nmod == 5){
                switchBit(n);
            }
        }
    }

    xroof = (unsigned) sqrt(((double)(limit - 1)) / 3);
    for(x = 1; x <= xroof; x++){
        yroof = (unsigned) sqrt((double)(limit - 3 * x * x));
        for(y = 1; y <= yroof; y++){
            n = (3 * x * x) + (y * y);
            nmod = n % 12;
            if (nmod == 7){
                switchBit(n);
            }
        }
    }

    xroof = (unsigned) sqrt(((double)(limit + 1)) / 3);
    for(x = 1; x <= xroof; x++){
        yroof = (unsigned) sqrt((double)(3 * x * x - 1));
        for(y = 1; y <= yroof; y++){
            n = (3 * x * x) - (y * y);
            nmod = n % 12;
            if (nmod == 11){
                switchBit(n);
            }
        }
    }
}


unsigned PrimeGen::next()
{
    while (bookmark <= limit)
    {
        bookmark++;

        if (getBit(bookmark))
        {
            unsigned out = bookmark;

            for(unsigned num = bookmark * 2; num <= limit; num += bookmark)
            {
                offBit(num);
            }

            return out;
        }
    }

    return 0;
}


inline void PrimeGen::printBoard()
{
    for(unsigned i = 0; i < sievelen; i++)
    {
        if (i % 4 == 0)
            cout << endl;

        printBinary(sieve[i]);
        cout << " ";
    }
}


inline unsigned char flag(const unsigned char bit_index)
{
    return ((unsigned char) 128) >> bit_index;
}


inline void printBinary(unsigned char byte)
{
    unsigned int i = 1 << (sizeof(byte) * 8 - 1);

    while (i > 0) {
        if (byte & i)
            cout << "1";
        else
            cout << "0";
        i >>= 1;
    }
}

我尽力清理它并使其可读。我不是专业程序员，所以请手下留情。

这是我得到的输出，当我初始化一个名为 pgen 的 PrimeGen 对象时，使用 pgen.printBoard() 打印其初始位板（请注意，在 next() 迭代之前缺少 23 和 59），然后迭代 next() 和打印所有返回的素数：

00000101 00010100 01010000 01000101
00000100 01010001 00000100 00000100
00010001 01000001 00010000 01000000
01000101 00010100 01000000 00000001

5
7
11
13
17
19
29
31
37
41
43
47
53
61
67
71
73
79
83
89
97
101
103
107
109
113
127

DONE

Process returned 0 (0x0)   execution time : 0.064 s
Press any key to continue.

Answer 1

尤里卡！！！

正如所料，这对我来说是一个愚蠢的错误。

3x^2 - y^2 方程有一个我忽略的小警告：x > y。考虑到这一点，我切换 23 和 59 的次数太多了，导致它们失败了。

感谢大家的帮助。救了我的培根。