素数求和问题,也是大一的一次实验。重新回顾,重新体会。
问题描述:从键盘输入任意一个整数n,编程计算并输出1~n之间所有素数之和。
附加题(选做):针对实验的问题想出一种算法,能对任意一个5<n<1,000,000,000的整数进行处理,并且时间限制在15秒内,内存利用限制为32M。
首先,必须了解下素数的概念: (百度百科) http://baike.baidu.com/view/10626.htm?fromId=1767
阶段一。常规逐个判断是否是素数
这里,n的大小没做具体要求,所用时间,所占内存也都没限制,所以代码比较随意,常规的判断一个数是否是素数,满足条件就累加的和上。
效率低,但也是我最初的方法。还是粘出来留念留念。
#include <stdio.h> #include <stdlib.h> #include <math.h> int is_prime( int n ) { int i, j, ret; int sum = 2; for ( j=2; j<=n; j++ ) { ret = 1; /* 利用ret的不同返回值,进行求和 */ if( j % 2 != 0 ) { for( i=3; i<=sqrt(j); i++ ) { if( j %i == 0 ) ret = 0; } } else ret = 0; if (ret == 1) sum = sum + j; } return sum; } int main() { int n, result; printf("please input a number:"); scanf("%d", &n); result = is_prime( n ); printf("%d", result); return 0; }
阶段二。素数筛选法。
素数筛选法之前我没接触过,没什么概念。是在尝试完成附加题的时候学习的。
基本思想
(以上内容来自百度百科)
现在,利用筛选法,写了个简单的例子,筛选出了100以内所有的素数。
#include <stdio.h> int main() { int a[101], i, j; for( i=2; i<=100; i++ ) //先初始化,且各位都不为0 a[i]=i; for(i=2;i<=50;i++) //直接判断 max/2即可 大于一半后该数的整数倍不在max范围内,无需考虑 { if(a[i]!=0) { for( j=i+i; j<=100; j+=i ) //每次把i的倍数置0,去掉。 a[j]=0; } } for(i=2;i<=100;i++) { if(a[i]!=0) //通过判断是否为0,进行筛选。 printf("%3d",a[i]); } return 0; }
相对来说,利用这个筛选法比常规的逐个判断在效率上有了很大的提高,可以说以我当时初学c的水平,最多也就只能写出这样的代码了。
可是,要想完成附加题,这样的显然是不行的。 且不说时间,就内存来说,32M = 33,554,432 < 1,000,000,000 * 4 。也不允许,所以,常规的筛法是不能完成的。
可是,附加题加分的阿。。 所以呢,,, 在万般无奈下,选择了百度.....就有了接下去的延伸。
阶段三。算法优化。
1.大牛算法一。 ----忘记出处了,在此道个歉了。
#include <time.h> #include <math.h> #define SIZE 1000000000 #define BITMAPSIZE 8 unsigned char map[SIZE / 2 / BITMAPSIZE + 1]; unsigned char bsm[10000];/* 这是一个估算值 */ int bn[10000]; int main() { int i, j, l; int c = 0; int st; int basn; st = clock(); /* 计算素数 */ basn = (int)sqrt(SIZE) + 1; for (i = 3; ; i += 2) { if ((bsm[i / BITMAPSIZE] & (1 << (i % BITMAPSIZE))) == 0) { for (j = i << 1; j < basn; j += i) { bsm[j / BITMAPSIZE] |= (1 << (j % BITMAPSIZE)); } bn[c++] = i; if (i > basn) break; } } for (i = 0; i < c; i++) { j = bn[ i ] << 1; l = j + bn[ i ]; do { map[l / (BITMAPSIZE * 2)] |= (1 << ((l >> 1) % BITMAPSIZE)); l += j; } while (l < SIZE); } c = 0; for (i = 1; i < SIZE / 2; i++) { i++; if ((map[i / BITMAPSIZE] & (1 << (i % BITMAPSIZE))) == 0) { c++; } i++; if ((map[i / BITMAPSIZE] & (1 << (i % BITMAPSIZE))) == 0) { c++; } } printf("1~1000000000 number: %d ", c); return 0; }
本质就是素数筛,只不过用了两遍。
第一组循环统计根号十亿内的素数。正常的筛法应用。
第二组循环用上面计算出来的素数筛剩余的范围。
这里比较有意思的地方就是它节省内存用的技巧。1.使用位作标志。2.它删掉了所有的偶数。也就是说第1位表示3,第2位表示5等等。也正因如此它的筛的步长才是两倍的素数,而起始步长是3倍的素数。
第三组循环就是统计标志位了。
2.大牛算法二。----bccn里面的beyondyf版主,杨大哥的算法。
在上一算法的基础上进一步优化,代码量是少了。同时,也更难懂了。
#include<stdio.h> #define RANGE 1000000000 char P[RANGE / 16 + 1]; int main() { int i, j, t, c = 1; for(i = 3; i <= RANGE; i += 2) if(!(P[i >> 4] & (1 << (i >> 1 & 7)))) for(c++, t = i + i, j = t + i; j > 0 && j <= RANGE; j += t) P[j >> 4] |= 1 << (j >> 1 & 7); printf("%d ", c); return 0; }
3.大牛算法三。------原文出处http://blog.csdn.net/redraiment/article/details/2072005 作者:子清行
算法优化
- 要找出一个数的因子,其实不需要检查 2→k,只要从 2->sqrt(k),就可以了。所有,我们筛法里,其实只要筛到sqrt(n)就已经找出所有的素数了,其中n为要搜索的范围。
- 另外,我们不难发现,每找到一个素数 k,就一次删除 2k, 3k, 4k,..., ik,不免还是有些浪费,因为2k已经在找到素数2的时候删除过了,3k已经在找到素数3的时候删除了。因此,当 i<k 时,都已经被前面的素数删除过了,只有那些最小的质因子是k的那些数还未被删除过,所有,就可以直接从 k*k 开始删除。
- 再有,所有的素数中,除了 2 以外,其他的都是奇数,那么,当 i 是奇数的时候,ik 就是奇数,此时 k*k+ik 就是个偶数,偶数已经被2删除了,所有我们就可以以2k为单位删除步长,依次删除 k*k, k*k+2k, k*k+4k, ...。
- 我们都清楚,在前面一小段范围内,素数是比较集中的,比如 1→100 之间就有25个素数。越到后面就越稀疏。
因为这些素数本身值比较小,所以搜索范围内,大部分数都是它们的倍数,比如搜索 1→100,这 100 个数。光是 2 的倍数就有 50 个,3 的倍数有 33 个,5的倍数 20 个,7 的倍数 14 个。我们只需搜索到7就可以,因此一共做删除操作50+33+20+14=117次,而 2 和 3 两个数就占了 83 次,这未免太浪费时间了。
所以我们考虑,能不能一开始就排除这些小素数的倍数,这里用 2 和 3 来做例子。
如果仅仅要排除 2 的倍数,数组里只保存奇数:1、3、5...,那数字 k 的坐标就是 k/2。
如果我们要同时排除 2 和 3 的倍数,因为 2 和 3 的最小公倍数是 6,把数字按 6 来分组:6n, 6n+1, 6n+2, 6n+3, 6n+4, 6n+5。其中 6n, 6n+2, 6n+4 是 2 的倍数,6n+3 是 3 的倍数。所以数组里将只剩下 6n+1 和 6n+5。n 从 0 开始,数组里的数字就一次是 1, 5, 7, 11, 13, 17...。
现在要解决的问题就是如何把数字 k 和它的坐标 i 对应起来。比如,给出数字 89,它在数组中的下标是多少呢?不难发现,其实上面的序列,每两个为一组,具有相同的基数 n,比如 1 和 5 ,同是 n=0 那组数,6*0+1 和 6*0+5;31 和 35 同是n=5那组,6*5+1 和 6*5+5。所以数字按6分组,每组2个数字,余数为5的数字在后,所以坐标需要加 1。
所以 89 在第 89/6=14 组,坐标为 14*2=28,又因为 89%6==5,所以在所求的坐标上加 1,即 28+1=29,最终得到 89 的坐标 i=29。同样,找到一个素数 k 后,也可以求出 k*k 的坐标等,就可以做筛法了。
这里,我们就需要用 k 做循环变量了,k 从 5 开始,交替与 2 和 4 相加,即先是 5+2=7,再是 7+4=11,然后又是 11+2=13...。这里我们可以再设一个变量gab,初始为 4,每次做 gab = 6 - gab,k += gab。让gab在2和4之间交替变化。另外,2 和 4 都是 2 的幂,二进制分别为10和100,6的二进制位110,所以可以用 k += gab ^= 6来代替。参考代码:
gab = 4; for (k = 5; k * k <= N; k += gab ^= 6) { ... }
但我们一般都采用下标 i 从 0→x 的策略,如果用 i 而不用 k,那应该怎么写呢?
由优化策略(1)可知,我们只要从 k2 开始筛选。 n=i/2,我们知道了 i 对应的数字 k 是素数后,根据(2),那如何求得 k2 的坐标 j 呢?这里假设 i 为偶数,即 k=6n+1。
k2 = (6n+1)*(6n+1) = 36n2 + 12n + 1,其中 36n2+12n = 6(6n2+2n) 是6的倍数,所以 k2 除 6 余 1。
所以 k2 的坐标 j = k2/6*2 = 12n2+4n。
由优化策略(2)可知,我们只要依次删除 k2+2l×k, l = 0, 1, 2...。即 (6n+1)×(6n+1+2l)。
我们发现,但l=1, 4, 7...时,(6n+1+2l)是3的倍数,不在序列中。所以我们只要依次删除 k2, k2+4l, k2+4l+2l...,又是依次替换2和4。
为了简便,我们可以一次就删除 k2 和 k2+4l 两项,然后步长增加6l。所以我们需要求 len=4l 和 stp=6l。不过这里要注意一点,k2+4k=(6n+1)*(6n+5),除以6的余数是5,坐标要加1。
len = k*(k+4)/6*2 - k
2
/6*2 = (6n+1)*(6n+1+4)/6*2+1 - (6n+1)*(6n+1)/6*2 = (12n
2
+12n+1) - (12n
2
+4n) = 8n+1; stp = k*(k+6)/6*2 - k
2
/6*2 = 12n+2;
最终,我们得到:
len = 8n+1; stp = 12n+2; j = 12n2+4n;
同理可以求出 k=6n+5 时的情况:
len = 4n+3; stp = 12n+10; j = 12n
2
+20n+8;
下面的代码在实现上用了位运算,可能有点晦涩。
#include <stdlib.h> #include <stdio.h> #include <time.h> #define N 1000000000 #define size (N/6*2 + (N%6 == 5? 2: (N%6>0))) int p[size / 32 + 1] = {1}; int creat_prime(void) { int i, j; int len, stp; int c = size + 1; for (i = 1; ((i&~1)<<1) * ((i&~1) + (i>>1) + 1) < size; i++) { if (p[i >> 5] >> (i & 31) & 1) continue; len = (i & 1)? ((i&~1)<<1) + 3: ((i&~1)<<2) + 1; stp = ((i&~1)<<1) + ((i&~1)<<2) + ((i & 1)? 10: 2); j = ((i&~1)<<1) * (((i&~1)>>1) + (i&~1) + 1) + ((i & 1)? ((i&~1)<<3) + 8 + len: len); for (; j < size; j += stp) { if (p[j >> 5] >> (j & 31) & 1 ^ 1) p[j >> 5] |= 1L << (j & 31), --c; if (p[(j-len) >> 5] >> ((j-len) & 31) & 1 ^ 1) p[(j-len) >> 5] |= 1L << ((j-len) & 31), --c; } if (j - len < size && (p[(j-len) >> 5] >> ((j-len) & 31) & 1 ^ 1)) p[(j-len) >> 5] |= 1L << ((j-len) & 31), --c; } return c; } int main(void) { clock_t t = clock(); printf("%d ", creat_prime()); printf("Time: %f ", 1.0 * (clock() - t) / CLOCKS_PER_SEC); return EXIT_SUCCESS; }
以上三种优化算法,都达到了要求。在我的机器上,算法3的耗时最短,仅仅4秒多。
不过就我个人目前的能力来说,看懂这些算法还是相当吃力的。 先做个标记,等以后回头了慢慢咀嚼。