素数筛选法优化_2.doc

2024年12月25日
素数筛选法优化
邓达生
原始的素数筛选法是在{2,……,n}里进行筛选。代码：
void workprime(unsigned int m, unsigned int *prime, unsigned int &nprime, char *mark){
unsigned int i, p;
nprime = -1;
memset(mark, 0, m);
for(i = 2; i < m; ++i)
if(!mark[i])
for(p = (prime[++nprime] = i) + i; p < m; p += i)
mark[p] = 1;
}

}

(注：本文章的代码都不进行for循环优化，例如上述代码改为
void workprime(unsigned int m, unsigned int *prime, unsigned int &nprime, char *mark){
unsigned int i, p, m2 = m >> 1;
nprime = -1;
memset(mark, 0, m);
for(i = 2; i <= m2; ++i)
if(!mark[i])
for(p = (prime[++nprime] = i) + i; p < m; p += i)
mark[p] = 1;
for(; i < m; ++i)
if(!mark[i])
prime[++nprime] = i;
}

若要优化请自行优化)
本文章的优化仅仅是加快筛选的速度(因为基础都是筛选)，如果原算法时间和空间复杂度是f(n)和g(n)，则优化后为k1f(n)和k2g(n)，k1、k2＜1。
优化：
因为2是已知的素数，所以{2,……,N}里的所有偶数(2的倍数)都可以在筛选前排除，也就是说待筛选列表里不应该存在2的倍数，并且在筛选中不应该存在检测2的倍数是否素数的行为或者将2的倍数标记为合数的行为，例如当检测到3是素数后，要把3的倍数都标记为合数时，不应该把6标记为合数，因为在待筛选列表里没有6。
我的基本思想就是在原始的待筛选列表里删除某(几)个已知素数的倍数。设删除后的待筛选列表为D，很容易证明若x、kx∈D，则k∈D。证明：若k∉D，则k是某个已知素数的倍数，则kx是某个已知素数的倍数，所以kx∉D，与已知kx∈D矛盾，所以k∈D。
现在原始待筛选列表删除2的倍数(即D={3,5,7,9……})，对于标记数组mark[]，mark[i]代表的就已经不能是i了，mark[i]代表的是2i+1是否合数。
筛选过程：假设现在mark[i]是false，代表2i+1是素数。然后就要把2i+1的倍数都标记为合数。应该被标记为合数的数分别是

应标记的数/2i+1的倍数
对应mark的下标
3(2i+1)=2(3i+1)+1
3i+1
5(2i+1)=2(5i+2)+1
5i+2
7(2i+1)=2(7i+3)+1
7i+3
……
……
(2k+1)(2i+1)=2((2k+1)i+k)+1
(2k+1)i+k

所以待标记的数的下标初始值为3i+1，增量为2i+1.
算法代码如下(m = n / 2)：void workprime(unsigned int m, unsigned int *prime, unsigned int &nprime, char *mark){
unsigned int i, j, p;
prime[nprime = 0] = 2;
memset(mark, 0, m);
for(i = 1; i < m; ++i)
if(!mark[i])
for(p = (prime[++nprime] = j = (i << 1) + 1) + i; p < m; p += j)
mark[p] = 1;
}

对于求n以内的素数，原始算法m=n，标记合数for循环增量是i，优化后的算法m=n/2，标记合数for循环增量是2i+1，效率提高大约一倍。
继续优化：原始待筛选列表删除2的倍数后m=n/2，再删除3的倍数后m=n*(2/6)=n/3，删除5的倍数后m=n*(8/30)=4n/15，……可以看出来，到后面即使再优化，效率也不会提高多少，而其后可以看到，效率在以很小的幅度提高时，代码的复杂度却以很高的速度提高。(以上为个人简单估计)
再删除3的倍数：D={5,7,11,13,17,19……}
令mark[i]代表6i±1是否合数。char类型是1字节8位，用低两位来标记(0位代表6i-1,1位代表6i+1)。

6i-1的倍数
mark下标
mark下标增量
位
6i-1
i

0
5(6i-1)=30i-5=6(5i-1)+1
5i-1
4i-1
1
7(6i-1)=42i-7=6(7i-1)-1
7i-1
2i
0
11(6i-1)=66i-11=6(11i-2)+1
11i-2
4i-1
1
13(6i-1)=78i-13=6(13i-2)-1
13i-2
2i
0
……
……
……
……
(6j-1)(6i-1)=6((6j-1)i-j)+1
(6j-1)i-j
4i-1
1
(6j+1)(6i-1)=6((6j+1)i-j)-1
(6j+1)i-j
2i
0

6i+1的倍数
mark下标
mark下标增量
位
6i+1
i

1
5(6i+1)=30i+5=6(5i+1)-1
5i+1
4i+1
0
7(6i+1)=42i+7=6(7i+1)+1
7i+1
2i
1
11(6i+1)=66i+11=6(11i+2)-1
11i+2
4i+1
0
13(6i+1)=78i+13=6(13i+2)+1
13i+2
2i
1
……
……
……
……
(6j-1)(6i+1)=6((6j-1)i+j)-1
(6j-1)i+j
4i+1
0
(6j+1)(6i+1)=6((6j+1)i+j)+1
(6j+1)i+j
2i
1
最终算法代码(m = n / 6)：void workprime(unsigned int m, unsigned int *prime, unsigned int &nprime, char *mark){
unsigned int i, j1, j2, p;
prime[0] = 2;
prime[nprime = 1] = 3;
memset(mark, 0, m);
for(i = 1; i < m; ++i){
if(!(mark[i] & 1)){
prime[++nprime] = (j1 = (i << 2) - 1) + (j2 = i << 1);
for(p = i + j1; p < m; p += j1){
mark[p] |= 2;
if((p += j2) >= m)break;
mark[p] |= 1;
}
}
if(!(mark[i] & 2)){
prime[++nprime] = (j1 = (i << 2) + 1) + (j2 = i << 1);
for(p = i + j1; p < m; p += j1){
mark[p] |= 1;
if((p += j2) >= m)break;
mark[p] |= 2;
};
};
};
};

再删除5的倍数：D={7,11,13,17,19,23,29,31……}，而mark[i]代表30i+1、30i+7、30i+11、30i+13、30i+17、30i+19、30i+23、30i+29 共8个数，刚刚一个字节可以存完，没有浪费任何空间。
注：除法 "/" 是整除，"%"运算是模运算(即求余)
30i+k的30j+c倍数
mark下标增量
(30j+1)(30i+k)=30((30j+1)i+jk)+k

(30j+7)(30i+k)=30((30j+7)i+jk+7k/30)+7k%30
6i+7k/30
(30j+11)(30i+k)=30((30j+11)i+jk+11k/30)+11k%30
4i+11k/30-7k/30
(30j+13)(30i+k)=30((30j+13)i+jk+13k/30)+13k%30
2i+13k/30-11k/30
(30j+17)(30i+k)=30((30j+17)i+jk+17k/30)+17k%30
4i+17k/30-13k/30
(30j+19)(30i+k)=30((30j+19)i+jk+19k/30)+19k%30
2i+19k/30-17k/30
(30j+23)(30i+k)=30((30j+23)i+jk+23k/30)+23k%30
4i+23k/30-19k/30
(30j+29)(30i+k)=30((30j+29)i+jk+29k/30)+29k%30
6i+29k/30-23k/30
(30(j+1)+1)(30i+k)=30((30j+31)i+jk+k)+k
2i+k-29k/30
这样的表一共有8个，共有8*8=64项，所以我之前说过代码复杂度会提高很多。如果就想前面那样写代码的话这个函数就会很长，所以就制表来解决。
代码很长，后面再列出来。为了测试上述3种算法的效率，我在中山大学的在线测评系统sicily的第一题"A-B"测试了运行效率，各自计算33554432(225)内的素数(共2063689个)，只删除2的倍数的用时大约0.9s，再删除3的倍数的用时大约0.6s，再删除5的倍数的用时只有0.2s！（全部代码在本机测试结果没错后再去sicily测试时间）难道是n不够大？还是我的估计误差太大呢？不过由于代码复杂度都已经比较高了，而且内存空间得到充分利用(mark的每一项代表8个数，刚好一个字节8位二进制用完)，我就不深入探讨了，有兴趣的读者可以自己算一下。下面附上代码(m = n / 30)：

void workprime(unsigned int m, unsigned int *prime, unsigned int &nprime, char *mark){
static const char marks[][8]={
2,4,8,16,32,64,128,1,
32,16,1,128,8,4,64,2,
16,1,64,2,128,8,32,4,
1,64,32,4,2,128,16,8,
128,2,4,32,64,1,8,16,
8,128,2,64,1,16,4,32,
4,8,128,1,16,32,2,64,
64,32,16,8,4,2,1,128};
static const unsigned int incs[][8]={
0,0,0,0,0,0,0,1,
1,1,1,0,1,1,1,1,
2,2,0,2,0,2,2,1,
3,1,1,2,1,1,3,1,
3,3,1,2,1,3,3,1,
4,2,2,2,2,2,4,1,
5,3,1,4,1,3,5,1,
6,4,2,4,2,4,6,1};
static const unsigned int base[]={1,7,11,13,17,19,23,29};
static const char basemark[]={1,2,4,8,16,32,64,128};
unsigned int i, i2, i4, i6, i30, p, j, k;
prime[0] = 2;
prime[1] = 3;
prime[2] = 5;
prime[3] = 7;
prime[4] = 11;
prime[5] = 13;
prime[6] = 17;
prime[7] = 19;
prime[8] = 23;
prime[nprime = 9] = 29;
memset(mark, 0, m);

for(k = 1; k < 8; ++k) // 标记 7、11、13、17、19、23、29 的倍数为合数(即 i = 0)
for(p = 0;;){
if((p += incs[k][0]) >= m)break; else mark[p] |= marks[k][0];
if((p += incs[k][1]) >= m)break; else mark[p] |= marks[k][1];
if((p += incs[k][2]) >= m)break; else mark[p] |= marks[k][2];
if((p += incs[k][3]) >= m)break; else mark[p] |= marks[k][3];
if((p += incs[k][4]) >= m)break; else mark[p] |= marks[k][4];
if((p += incs[k][5]) >= m)break; else mark[p] |= marks[k][5];
if((p += incs[k][6]) >= m)break; else mark[p] |= marks[k][6];
if((p += incs[k][7]) >= m)break; else mark[p] |= marks[k][7];
};
for(i = 1; i < m; ++i){
i6 = (i2 = i << 1) + (i4 = i << 2);
i30 = (i << 5) - i2;
for(k = 0; k < 8; ++k)
if(!(mark[i] & basemark[k])){
prime[++nprime] = i30 + base[k];
for(p = i;;){
if((p += i6 + incs[k][0]) >= m)break; else mark[p] |= marks[k][0];
if((p += i4 + incs[k][1]) >= m)break; else mark[p] |= marks[k][1];
if((p += i2 + incs[k][2]) >= m)break; else mark[p] |= marks[k][2];
if((p += i4 + incs[k][3]) >= m)break; else mark[p] |= marks[k][3];
if((p += i2 + incs[k][4]) >= m)break; else mark[p] |= marks[k][4];
if((p += i4 + incs[k][5]) >= m)break; else mark[p] |= marks[k][5];
if((p += i6 + incs[k][6]) >= m)break; else mark[p] |= marks[k][6];
if((p += i2 + incs[k][7]) >= m)break; else mark[p] |= marks[k][7];
};
};
};
};

对于标记i的倍数ji为合数时，只需要从j=i，即i2开始标记，因为对于j<i，ji一定已经标记了，当j的质因数被找到。所以前面的代码可以再优化一下，因为i2很容易超过数据范围所以要先√n一下。下面代码到sicily测试(前面的一样)，用时0.13s。
unsigned int workprime(unsigned int m, unsigned int *prime, char *mark){
unsigned int nprime;
static const char marks[][8]={
2,4,8,16,32,64,128,1,
32,16,1,128,8,4,64,2,
16,1,64,2,128,8,32,4,
1,64,32,4,2,128,16,8,
128,2,4,32,64,1,8,16,
8,128,2,64,1,16,4,32,
4,8,128,1,16,32,2,64,
64,32,16,8,4,2,1,128
};
static const unsigned int incs[][8]={
0,0,0,0,0,0,0,1,
1,1,1,0,1,1,1,1,
2,2,0,2,0,2,2,1,
3,1,1,2,1,1,3,1,
3,3,1,2,1,3,3,1,
4,2,2,2,2,2,4,1,
5,3,1,4,1,3,5,1,
6,4,2,4,2,4,6,1
};
static const unsigned int base[]={1,7,11,13,17,19,23,29};
static const char basemark[]={1,2,4,8,16,32,64,128};
unsigned int i, i2, i4, i6, i30, p, k, m2 = 1 + (int)sqrt(m / 30.0);
prime[0] = 2;
prime[1] = 3;
prime[2] = 5;
prime[3] = 7;
prime[4] = 11;
prime[5] = 13;
prime[6] = 17;
prime[7] = 19;
prime[8] = 23;
prime[nprime = 9] = 29;
memset(mark, 0, m);

for(k = 1; k < 8; ++k)
for(p = 0;;){
if((p += incs[k][0]) >= m)break; else mark[p] |= marks[k][0];
if((p += incs[k][1]) >= m)break; else mark[p] |= marks[k][1];
if((p += incs[k][2]) >= m)break; else mark[p] |= marks[k][2];
if((p += incs[k][3]) >= m)break; else mark[p] |= marks[k][3];
if((p += incs[k][4]) >= m)break; else mark[p] |= marks[k][4];
if((p += incs[k][5]) >= m)break; else mark[p] |= marks[k][5];
if((p += incs[k][6]) >= m)break; else mark[p] |= marks[k][6];
if((p += incs[k][7]) >= m)break; else mark[p] |= marks[k][7];
};
for(i = 1; i <= m2; ++i){
i6 = (i2 = i << 1) + (i4 = i << 2);
i30 = (i << 5) - i2;
for(k = 0; k < 8; ++k)
if(!(mark[i] & basemark[k])){
prime[++nprime] = i30 + base[k];
if((p = (i30 + 1) * i + base[k] * i) >= m)continue;
for(mark[p] |= marks[k][7];;){
if((p += i6 + incs[k][0]) >= m)break; else mark[p] |= marks[k][0];
if((p += i4 + incs[k][1]) >= m)break; else mark[p] |= marks[k][1];
if((p += i2 + incs[k][2]) >= m)break; else mark[p] |= marks[k][2];
if((p += i4 + incs[k][3]) >= m)break; else mark[p] |= marks[k][3];
if((p += i2 + incs[k][4]) >= m)break; else mark[p] |= marks[k][4];
if((p += i4 + incs[k][5]) >= m)break; else mark[p] |= marks[k][5];
if((p += i6 + incs[k][6]) >= m)break; else mark[p] |= marks[k][6];
if((p += i2 + incs[k][7]) >= m)break; else mark[p] |= marks[k][7];
};
};
};
for(; i < m; ++i){
i6 = (i2 = i << 1) + (i4 = i << 2);
i30 = (i << 5) - i2;
for(k = 0; k < 8; ++k)
if(!(mark[i] & basemark[k]))
prime[++nprime] = i30 + base[k];
};
return nprime;
};

(对这篇文章有什么建议或问题可发邮件到，谢谢)
10