#include<stdio.h>
#include<time.h>
#include<stdlib.h>
#define m 10000001
int main()
{
clock_t start, finish;
int i,j,p[m],count;
start=clock();
for(i=0;i<m;i++)
p[i]=1;
for(i=2;i*i<m;i++)
if(p[i])
{
for(j=i+i;j<m;j+=i)
p[j]=0;
}
count=0;
for(i=2;i<m;i++)
if(p[i])
count++;
finish=clock();
printf("%f.3",(double)(finish - start) / CLOCKS_PER_SEC);
return 0;
}
质数算法大全,及C程序实现优化详解 筛选法
三木追风
质数的定义
一个数,如果只有1和它本身两个因数,这样的数叫做质数,又称素数。
在上文 《素数算法大全,及C程序实现优化详解 (一) 试除法》中我们已经探讨了求解素数的一类算法,并且将试除法从最初的低效版本优化的高效的V2。那么,还有没有其它更佳算法呢?这就是下面三藏要和大家探讨的内容
合数过滤筛选法
算法描述:我们知道,素数N不能被2~(N-1)间的任何数整除;反过来看,只要能被2~(N-1)间的任何数整除的N,都不是素数。所以我们可以采用一个简单的排除法:就是对N以内的所有数,只要逐个去除值为2~(N-1)的倍数的数,剩下的就是素数。
C语言实现
// 合数过滤筛选法 Ver1
// 参数:n 求解n以内(包括n)的素数
// 返回值:n以内素数个数
int CompositeNumFilterV1(int n)
{
int i, j;
// 素数数量统计
int count = 0;
// 分配素数标记空间,结合后文思考为何+1
char* flag = (char*)malloc( n+1 );
// 初始化素数标记
for (i=2; i<=n; i++)
{
// 为什么*(p+i)要写成flag[i]呢?可读性更佳尔
flag[i] = 1;
}
// 写程序要注意排版和留空,方便阅读,也可减少出错几率
// 以2~(N-1)为因子过滤合数
for (i=2; i < n; i++)
{
for (j=2; i*j <= n; j++)
{
// i*j是由i,j两整数相乘而得,显然不是素数
flag[i*j] = 0;
}
}
// 统计素数个数
for (i=2; i<=n; i++)
{
// 其实if(flag)就其同样作用了,但这么写是有留言的
// 请参阅《C语言程序设计常见错误剖析及解决之道》一文
if (1 == flag[i]) count++;
}
// 因输出费时,且和算法核心相关不大,故略
// 释放内存,别忘了传说中的内存泄漏
free(flag);
return count;
}
在上文给出的main函数中以不同参数调用CompositeNumFilterV1函数,得到执行结果如下:
[100000]以内素数个数:9592, 计算用时:15毫秒
[1000000]以内素数个数:78498, 计算用时:125毫秒
[5000000]以内素数个数:348513, 计算用时:2578毫秒
[10000000]以内素数个数:664579, 计算用时:6281毫秒
注:因程序是非独占性运行的,所以时间不是完全精确的,但基本能反映实情
显然,比上文中的试除法要快,而且谁都可以看到上例是一个未经优化的粗陋版本,好多地方是三藏故意采用比较低效做法,为了与后文的优化版比较,凸显优化之重要,也为了初学者记住别采用类似低效做法,下面我们开始优化之旅
优化分析
上面CompositeNumFilterV1函数存在的问题有:
1. 在外层循环,需要一直执行到n-1吗?不要,因为n/2~n-1间的数显然不能整出n
2. 在内层循环中重复使用i*j显然是低效的,考虑到计算机中加减运算速度比乘除快,可以考虑变乘法为加法
3. 在循环修改flag过程中,其实有很多数会被重复计算若干次,比如6=2*3=3*2,会被重复置0,类似操作很多,所以我们得设法避免或减少flag重复置0
据上述分析,我们可将程序优化如下:
// 合数过滤筛选法 Ver2
// 参数:n 求解n以内(包括n)的素数
// 返回值:n以内素数个数
int CompositeNumFilterV2(int n)
{
int i, j;
// 素数数量统计
int count = 0;
// 分配素数标记空间,明白+1原因了吧,因为浪费了一个flag[0]
char* flag = (char*)malloc( n+1 );
// 初始化素数标记,要高效点咯
flag[2] = 1;
// 注意是i<n不是上例中的i<=n了,理由自思
for (i=3; i<n; i++)
{
flag[i++] = 1;
// 偶数自然不是素数,直接置0好了
flag[i] = 0;
}
// n为奇数
if (n%2 != 0)
{
flag[n] = 1;
}
// 从3开始filter,因为2的倍数早在初始化时代就干掉了
// 到n/2止的理由还要说吗
for (i=3; i <= n/2; i++)
{
// i是合数,请歇着吧,因为您的工作早有您的质因子代劳了
if (0 == flag[i]) continue;
// 从i的2倍开始过滤,变乘法为加法
for (j=i+i; j <= n; j+=i)
{
flag[j] = 0;
}
}
// 统计素数个数
for (i=2; i<=n; i++)
{
if (flag[i]) count++;
}
// 因输出费时,且和算法核心相关不大,故略
// 释放内存,别忘了传说中的内存泄漏
free(flag);
return count;
}
再来调用CompositeNumFilterV2得到执行结果:
[100000]以内素数个数:9592, 计算用时:n太小,时间精度不够
[1000000]以内素数个数:78498, 计算用时:31毫秒
[5000000]以内素数个数:348513, 计算用时:453毫秒
[10000000]以内素数个数:664579, 计算用时:1062毫秒
[100000000]以内素数个数:5761455, 计算用时:12973毫秒
哇哇,比昨天的试除发快了好多倍,可见算法的威力,值得好好学习,别说学算法没用咯。
上例着那个计算一亿以内的素数只要约13秒,应该算不错了,今天是否可以休息了呢?No,我们要追求极限!
int CompositeNumFilterV3(int n)
{
int i, j;
// 素数数量统计
int count = 0;
// 分配素数标记空间,明白+1原因了吧,因为浪费了一个flag[0]
char* flag = (char*)malloc( n+1 );
// 干嘛用的,请仔细研究下文
int mpLen = 2*3*5*7*11*13;
char magicPattern[mpLen];
// 奇怪的代码,why,思考无法代劳,想!
for (i=0; i<mpLen; i++)
{
magicPattern[i++] = 1;
magicPattern[i++] = 0;
magicPattern[i++] = 0;
magicPattern[i++] = 0;
magicPattern[i++] = 1;
magicPattern[i] = 0;
}
for (i=4; i<=mpLen; i+=5)
magicPattern[i] = 0;
for (i=6; i<=mpLen; i+=7)
magicPattern[i] = 0;
for (i=10; i<=mpLen; i+=11)
magicPattern[i] = 0;
for (i=12; i<=mpLen; i+=13)
magicPattern[i] = 0;
// 新的初始化方法,将2,3,5,7,11,13的倍数全干掉
// 而且采用memcpy以mpLen长的magicPattern来批量处理
int remainder = n%mpLen;
char* p = flag+1;
char* pstop = p+n-remainder;
while (p < pstop)
{
memcpy(p, magicPattern, mpLen);
p += mpLen;
}
if (remainder > 0)
{
memcpy(p, magicPattern, remainder);
}
flag[2] = 1;
flag[3] = 1;
flag[5] = 1;
flag[7] = 1;
flag[11] = 1;
flag[13] = 1;
// 从17开始filter,因为2,3,5,7,11,13的倍数早被kill了
// 到n/13止的,哈哈,少了好多吧
int stop = n/13;
for (i=17; i <= stop; i++)
{
// i是合数,请歇着吧,因为您的工作早有您的质因子代劳了
if (0 == flag[i]) continue;
// 从i的17倍开始过滤
int step = i*2;
for (j=i*17; j <= n; j+=step)
{
flag[j] = 0;
}
}
// 统计素数个数
for (i=2; i<=n; i++)
{
if (flag[i]) count++;
}
// 因输出费时,且和算法核心相关不大,故略
// 释放内存,别忘了传说中的内存泄漏
free(flag);
return count;
}
再看CompositeNumFilterV3执行结果:
[1000000]以内素数个数:78498, 计算用时:15毫秒
[5000000]以内素数个数:348513, 计算用时:203毫秒
[10000000]以内素数个数:664579, 计算用时:515毫秒
[100000000]以内素数个数:5761455, 计算用时:6421毫秒
再次优化后速度提升了又一倍左右,三藏不禁有点满足了,睡觉也!
素数筛选法优化
邓达生
原始的素数筛选法是在{2,……,n}里进行筛选。代码:
(注:本文章的代码都不进行for循环优化,例如上述代码改为
若要优化请自行优化)
本文章的优化仅仅是加快筛选的速度(因为基础都是筛选),如果原算法时间和空间复杂度是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):
对于求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):
再删除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):