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

注:这是本人去年非典期间在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后缀