素数算法逐步优化

柔光的暖阳◎ 2021-11-22 16:31 382阅读 0赞

素数求和问题，也是大一的一次实验。重新回顾，重新体会。

**问题描述：从键盘输入任意一个整数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;
    }

# 阶段二。素数筛选法。 #

素数筛选法之前我没接触过，没什么概念。是在尝试完成附加题的时候学习的。

### 基本思想 ###

用筛法求素数的基本思想是：把从1开始的、某一范围内的正整数从小到大顺序排列， 1不是素数，首先把它筛掉。剩下的数中选择最小的数是素数，然后去掉它的倍数。依次类推，直到筛子为空时结束。如有：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

1不是素数，去掉。剩下的数中2最小，是素数，去掉2的倍数，余下的数是：

3 5 7 9 11 13 15 17 19 21 23 25 27 29

剩下的数中3最小，是素数，去掉3的倍数，如此下去直到所有的数都被筛完，求出的素数为：

2 3 5 7 11 13 17 19 23 29

（以上内容来自百度百科）

现在，利用筛选法，写了个简单的例子，筛选出了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\n", 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\n", c);
        return 0;
    }

## 3.大牛算法三。------原文出处http://blog.csdn.net/redraiment/article/details/2072005 作者：子清行 ##

### 算法优化 ###

1.  要找出一个数的因子，其实不需要检查 2→k，只要从 2->sqrt(k)，就可以了。所有，我们筛法里，其实只要筛到sqrt(n)就已经找出所有的素数了，其中n为要搜索的范围。
2.  另外，我们不难发现，每找到一个素数 k，就一次删除 2k, 3k, 4k,..., ik，不免还是有些浪费，因为2k已经在找到素数2的时候删除过了，3k已经在找到素数3的时候删除了。因此，当 i＜k 时，都已经被前面的素数删除过了，只有那些最小的质因子是k的那些数还未被删除过，所有，就可以直接从 k\*k 开始删除。
3.  再有，所有的素数中，除了 2 以外，其他的都是奇数，那么，当 i 是奇数的时候，ik 就是奇数，此时 k\*k+ik 就是个偶数，偶数已经被2删除了，所有我们就可以以2k为单位删除步长，依次删除 k\*k, k\*k+2k, k\*k+4k, ...。
4.  我们都清楚，在前面一小段范围内，素数是比较集中的，比如 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)  
    {  
        ...  
    }

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 - k2/6*2 = (6n+1)*(6n+1+4)/6*2+1 - (6n+1)*(6n+1)/6*2 = (12n2+12n+1) - (12n2+4n) = 8n+1;
    stp = k*(k+6)/6*2 - k2/6*2 = 12n+2;

最终，我们得到：

len = 8n+1;
    stp = 12n+2;
      j = 12n2+4n;

同理可以求出 k=6n+5 时的情况：

len = 4n+3;
    stp = 12n+10;
      j = 12n2+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 \n", creat_prime());
        printf("Time: %f ", 1.0 * (clock() - t) / CLOCKS_PER_SEC);
        return EXIT_SUCCESS;
    }

以上三种优化算法，都达到了要求。在我的机器上，算法3的耗时最短，仅仅4秒多。

不过就我个人目前的能力来说，看懂这些算法还是相当吃力的。 先做个标记，等以后回头了慢慢咀嚼。