HDOJ 1164 Eddy's research I
质数的定义
一个数，如果只有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毫秒
再次优化后速度提升了又一倍左右，三藏不禁有点满足了，睡觉也！
1 #include<iostream> 2 //#include<algorithm> 3 const int MAX= 65536; 4 using namespace std; 5 6 int prime[MAX+1]; 7 int flag[MAX]; 8 9 void Get_Prime(int &count)//得到0~MAX内的所有素数 10 { 11 int i,j; 12 //memset(flag,true,sizeof(flag)); 13 flag[2]=1; 14 for(i=3;i<MAX;i++) 15 { 16 flag[i++]=true; 17 flag[i]=false;//所有偶数都不是素数 18 } 19 // MAX为奇数 20 if (MAX%2 != 0)//MAX不是偶数 21 { 22 flag[MAX] = true; 23 } 24 for(i=3 ;i<MAX/2; i++)//在外层循环，不需要一直执行到MAX 25 //因为MAX/2~MAX-1间的数显然不能整除MAX 26 { 27 if(flag[i]== false) continue; 28 for(j=2;i*j<MAX;i++) 29 { 30 flag[i*j]==false; 31 } 32 } 33 34 for(i=2;i<MAX;i++) 35 { 36 if(flag[i]) 37 { 38 prime[count]=i; 39 count++; 40 } 41 } 42 43 } 44 int main() 45 { 46 int cases; 47 int count=0; 48 Get_Prime(count); 49 while(scanf("%d",&cases)!=EOF) 50 { 51 int flag=true; 52 for(int i=0;i<count;i++) 53 { 54 while(cases%prime[i]==0) 55 { 56 if(flag) 57 { 58 flag=false; 59 cout<<prime[i]; 60 } 61 else 62 cout<<"*"<<prime[i]; 63 cases/=prime[i]; 64 } 65 if(cases==1) break; 66 } 67 cout<<endl; 68 } 69 70 }
