CSDN博客

img l1t

用C语言思想改写的用筛法求质数程序(第2修订版)

发表于2004/3/2 10:25:00  4363人阅读

注:这是本人去年非典期间在vccode上发表的拙作,近日看到这里人气不错,重贴给同好探讨,勿怪

http://www.vccode.com/file_show.php?id=1852

http://www.vccode.com/file_download.php?id=1852 源代码

原创: 本文及程序可免费使用,请勿用于商业用途

原算法思想:
1)先假定所有数都是质数,直到证明它不是。
2)通过筛除所有比较小的质数的倍数来找到质数。

改进一:
利用“大于2的质数都是奇数”这一知识,将存放待筛整数空间减少一半,只存储奇数。

改进二:
由于不再存有偶数,筛出的必然是奇数和奇数之积,乘数从3开始每次增2,将循环次数再减少一半

改进三:
(可能只能在C语言里这么做),利用位(bit)存放各数是否质数的标志,将存放空间减少到1/8
虽然每次循环的位运算步骤增加了,但是申请内存的时间大大减少,而且能求出更大的质数
定义了2个宏
#define getisp(i) (isprime[(i)/8]>>((i)%8))&0x1
用来取得第i个数是否质数标记
#define setisp(i,j) j?
(isprime[(i)/8]|=(0x1<<((i)%8))):
(isprime[(i)/8]&=~(0x1<<((i)%8)))
用来设置第i个数是否质数标记
不用函数表示的原因是为了较少的堆栈操作

改进四:
在设置第i个数的标记前先判断它是否已经设置为合数,如果已经设置为合数,则不重新设置
由于读操作快于写操作,这个改进对速度提高关系很大

改进五:
(1)把位运算宏中的/8操作改为>>3,把%8操作改为&0x7,
(2)由于只有一处是把数3的标记设为1,而且也不必要,把setisp宏分解,省略了设置值是1或0的判断,
(3)为了减少除法操作次数,把循环变量的值作了代数变换,但程序的可读性下降很多
这些改进对速度提高关系不大

改进六:
(1)每次筛除某个质数的倍数时,因为乘数小于该质数的积都已经被作为小于该质数的积被筛除了,
因此倍数应该从该质数开始,如:
筛除5的倍数,15已经作为3的倍数被筛除了,应该从25开始。这是算法的改进。
(2)考虑到位操作取得第i个数是否质数标记由4步运算组成,而用一个中间变量判断3的倍数只需要
3步运算。(一次自增,一次判断,判断成立再赋值一次),把判断3的倍数单独提出来。
至于大于3的质数的倍数,由于单独判断的语句多于用宏判断,得不偿失,故不再额外处理。
上述算法中其实3的倍数的存储空间也是可以节省的,但处理语句复杂化,速度上得不偿失,故也不再额外处理。

替换的核心语句如下:
由于程序的可读性原因附上了功能相同的语句作为注释:

      char _3times=0;
      i=(max-1)/2;
      for(int j = 4; j <=i ; j +=3) //set all 3’s times
      {
         isprime[(j)>>3]&=~(0x1<<((j)&0x7));
      }
      char i1=0;
       int m= (n-1)/2;
      int k=(max-1)/2;
      for( i = 2; i <= m; i++)
      {
         if(++i1==3) //skip when i is 3’s times
         {
           i1=0;
           continue;
         }
        if(getisp(i))
        {
           _3times=((i<<2)+1)%3;  //(2*(2*i+1)%3-1)

           for(int j = i*(2*i+2); j <= k; j = j +i +i+1)
            {
              if(++_3times==3) //skip when j is 3’s times
             {
                // printf("skip i=%d,j=%d%,%d
",2*i+1,2*j+1,_3times);
                _3times=0;
                continue;
             }
              if(getisp(j)) //good idea
               isprime[(j)>>3]&=~(0x1<<((j)&0x7));
             //else
              //ct++;
            }
          }
       }
      /*
       for( i = 5; i <= n; i+=2)
       {
         if(++i1==3) //skip when i is 3’s times
         {
           i1=0;
           continue;
         }
         int m=(i-1)/2;

         if (getisp(m))      // If i is a prime,
         {
             _3times=(i+i)%3-1; //用于定位j第几次是3的倍数(nod * nod )%3 all ==1
             for(int j = i*i; j <= max; j = j +i +i) // loop through multiples,
             {
               if(++_3times==3) //skip when j is 3’s times
               {
                 // printf("skip i=%d,j=%d%,%d
",i,j,_3times);
                 _3times=0;
                 continue;
               }
               int k=(j-1)/2;// they are not prime.
               if(getisp(k)) //good idea
               {
                     isprime[(k)>>3]&=~(0x1<<((k)&0x7)); // they are not prime.
               }
               else
               {
                 //printf("cat i=%d,j=%d%,%d
",i,j,_3times);
                 ct++;
               }
             }
           }
         }
       */
改进七:
昨天读了Java语法,发现它也支持位运算,但不支持宏定义,有更清晰的类型检查,如整数不能默认转化为boolean
于是,在改进六C程序的基础上,把isprime的类型 unsigned char* 改为 byte[],同时作为类数据成员变量。

static byte[] isprime;
把getisp改成如下的函数:
public static int getisp(int i)
    {
    return (isprime[(i)>>3]>>((i)&0x7))&0x1;
    }

调用时用if(getisp(j)==1)代替if(getisp(j)).
虽然Java的byte类型是有符号的,但没有影响程序的正确运行。
结果令人振奋,运行速度达到了改进3的C程序的水平。考虑到函数调用的开销,如果都改为展开的语句,
由于Java不支持if((isprime[(k)>>3]>>((k)&0x7))&0x1)==1)这样的表达式,必须先把位运算的结果赋值给某int变量f,
再判断表达式f==1是否成立,结果速度变化不明显。
补充:
利用ms java sdk 4.0提供的jexegen把.class文件转化为exe文件(但要求.class必须使用Ms的jvc而不能用Sun Jdk的javac编译,而且只能用到
jdk 1.1的功能)可提高Java程序速度,不过这已经不是纯粹的Java程序了
C:sieve>jexegen  /out:Sieve4.exe /main:Sieve4 Sieve4.class

测试速度方法:
命令行方式,用批处理test 程序名 最大整数的格式可看到程序执行前后的时间
例如:test java Sieve 50000000
注:这个时间包括class调入内存,解释的时间,对Java不太公平
测试结果:
平台:PIII650 128M windows 2003 jdk 1.4.1_02 lccwin32 ver3.3
C程序测试命令行
C:>test SieveX 12345678
Java程序测试命令行
C:>test Java SieveX 12345678
运算结果都是:
The largest prime less than or equal to 12345678 is 12345653,
平均用时如下:
/*------------原始java--------7.08 s------*/

/*------------改进一、二java------3.02 s--------*/

/*------------改进一、二c------2.34 s--------*/

/*------------改进三--------1.51 s------*/

/*------------改进四、五------1.12 s--------*/

/*------------改进六java------2.30 s--------*/
/*------------改进六java exe------1.64 s--------*/

/*------------改进七java------1.53 s--------*/
/*------------改进七java exe------1.14 s--------*/
结论:
1.算法的改进远比语句或指令的改进对提高程序的运行效率重要,除法指令没有我们想象中慢;
2.Java和C程序的运行速度有时候很接近;

C语言我能做到的改进就这么些了,下一步该用汇编了,可是我不会写汇编程序:(
数学上有没有更好的求质数的定理可用呢?我不知道,欢迎大家指正。
因为本程序主要用来说明求质数的问题,有些细节没有考虑,如最大数是0、1时的处理。
还有些地方,如x*3是否改成(x<<1)+x,j=j+i+i+1是否改成j=j+1+(i<<1),我觉得
不了解编译器的优化情况,就不自作多情了。
说明:
由于受int类型的取值范围和malloc空间限制,程序1-2能计算的质数不应超过10^8
程序3-5能计算的质数不应超过2*10^8,否则可能造成系统崩溃。

附件:
原Java程序 Sieve.java
你在编译时可以去掉package com.davidflanagan.examples.basics;一行,否则你需要设置
CLASSPATH环境变量,并把java源程序放在特定的目录才能用java com.davidflanagan.examples.basics.Sieve 2
去掉上述行的Java程序 Sieve1.java
包括改进一、二的Java程序 Sieve2.java
从改进Java程序改写的C程序 Sieve2.c
包括改进三的C程序 Sieve3.c
包括改进四、五的C程序 Sieve4.c
包括改进四、五、六的C程序 Sieve5.c
包括改进一、二、六(1)的java程序 Sieve4.java
包括改进七的java程序 Sieve5.java
批处理 test.bat
注:为了在VC下通过编译,需要用cpp后缀来表示C程序,在lcc下需要用c后缀

0 0

相关博文

我的热门文章

img
取 消
img