最近我一直在研究使用Atkin筛法（http://en.wikipedia.org/wiki/Sieve_of_atkin）生成质数的C++程序。我的目标是能够生成任何32位数字，主要用于欧拉计划问题。这只是一个暑期项目。
该程序使用位板来存储素性：即一系列由1和0组成的数字，例如第11位为1，第12位为0，第13位为1等。为了有效地使用内存，实际上这是一个char数组，每个char包含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()（请注意在下一次迭代之前缺少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.