zoukankan      html  css  js  c++  java
  • 素数算法逐步优化

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


    问题描述:从键盘输入任意一个整数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
    ", 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    作者:子清行

    算法优化

    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)  
    {  
        ...  
    }  

    但我们一般都采用下标 i 从 0→x 的策略,如果用 i 而不用 k,那应该怎么写呢?

    由优化策略(1)可知,我们只要从 k2 开始筛选。 n=i/2,我们知道了 i 对应的数字 k 是素数后,根据(2),那如何求得 k的坐标 j 呢?这里假设 i 为偶数,即 k=6n+1。

    k2 = (6n+1)*(6n+1) = 36n2 + 12n + 1,其中 36n2+12n = 6(6n2+2n) 是6的倍数,所以 k除 6 余 1。

    所以 k的坐标 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。

    为了简便,我们可以一次就删除 k和 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秒多。

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



  • 相关阅读:
    2020.10.23 19级training 补题报告
    2020.10.17 天梯赛练习 补题报告
    2020.10.16 19级training 补题报告
    2020.10.9 19级training 补题报告
    2020.10.10 天梯赛练习 补题报告
    2020.10.3 天梯赛练习 补题报告
    2020.10.2 19级training 补题报告
    第十届山东省ACM省赛复现补题报告
    VVDI Key Tool Plus Adds VW Passat 2015 Key via OBD
    Xhorse VVDI Prog Software V5.0.3 Adds Many MCUs
  • 原文地址:https://www.cnblogs.com/riskyer/p/3233899.html
Copyright © 2011-2022 走看看