筛法“四不像”——《C解毒》试读

题目：用筛法求100之内的素数。

#include <stdio.h>
#include <math.h>//程序中用到求平方根函数sqrt 
int main( )
{int i,j,n,a[101];                //定义a数组包含101个元素 
for(i=1;i<=100;i++)             //a[0]不用，只用a[1]~a[100] 
      a[ i ]=i;                     //使a[1]~a[100]的值为1~100 
  a[1]=0;                         //先“挖掉”a[1]
for(i=2;i<sqrt(100);i++)
     for(j=i+1;j<=100;j++)
        {if(a[ i ]!=0&&a[j]!=0)
            if(a[j]%a[ i ]==0)
              a[j]=0;             //把非素数“挖掉”
        } 
  printf("\n" )  ;
  for(i=2,n=0;i<=100;i++)
    { if(a[ i ]!=0)                 //选出值不为0的数组元素，即素数
        {printf("%5d",a[ i ]);      //输出素数，宽度为5列 
         n++;                     //累计本行已输出的数据个数 
        }
      if(n==10)
        {printf("\n" ) ;
         n=0;
        }
    }    
  printf("\n" ) ;  
  return0;
}

——谭浩强 ，《C程序设计（第四版）学习辅导》，清华大学出版社，2010年7月，p54~55

    这段代码的基本出发点是，用

int a[101];

这个数组来存储1~100这100个正整数，然而a却是一个有101个元素的数组。这种定义时空出一个多余的存储空间的做法非但没有任何必要性，而且有着双重的不必要性。
    首先，压根就没必要存储1~100这100个整数。但为了说明该程序中存在的其他问题，稍后的讨论暂时还是在存储1~100这100个整数的基础上进行。至于为什么不需要存储这100个数留在把这段代码中的其他问题讨论完毕之后再讨论。
    其次，即使采用存储1~100这100个正整数的方案来写代码，那么从1到100一共有100个整数，存储这100个整数也根本并不需要定义一个拥有101个int类型元素的数组（参见§55，p33），而应该定义为:

int a[100];

    这样，对该数组赋初值部分应写为

for ( i =0 ; i <100 ; i ++ )
     a [ i ] = i +1;

    亦即把1放在a[0]中，把2放在a[1]中……把100放在a[99]中。
    可能有人会辩解说定义有101个元素的数组可以把1放在a[1]中，把2放在a[2]中……把100放在a[100]中，代码似乎更容易写。然而这种辩解既是苍白无力也是毫无道理的。试想，如果问题要求存储10001~20000这10000个数，难道就会为此而定义一个尺寸为20001大小的数组，然后只使用其中的一半而把另一半不管不顾地浪费掉吗？
    下面这行代码

for(i=2;i<sqrt(100);i++)

存在的问题更多：
1.        首先，表达式i<sqrt(100)中的“<”其实应该是“<=”。为了说明这一点，首先来看一下究竟什么是筛法。
    筛法是公元前250年由古希腊数学家埃拉托斯特尼所提出的一种求出不超过自然数n范围内所有素数的一种简单的古老算法。具体做法是，首先写出1~n这n个自然数，这里假设n等于9：
1 2 3 4 5 6 7 8 9
由于1不是素数，所以把它划去
1 2 3 4 5 6 7 8 9
这样2就是找到的第一个素数，再把2的倍数划去
1 2 3 4 5 6 7 8 9
2之后第一个没被划去的3就是第二个素数，再把3的倍数划去
1 2 3 4 5 6 7 8 9
由于3等于9的平方根，所以到此时为止没被划去的数就构成了不超过自然数n范围内的所有素数。因为任何一个合数n，都必定存在一个小于或等于（而非小于）√n 且不等于1的因子。所以只要用所有的小于等于√n的素数去“筛”这张表得到的一定是1~n范围的所有素数。但是，如果只用所有小于√n的素数去“筛”这张表则是不充分的。例如，若前面只把2的倍数筛去，剩下的数就不全是素数。
    把题目中的100换成121，101换成122，亦即求1~121之间的所有素数，立刻就会发现输出结果是错误的，而这个错误的原因就在于代码中错误地把“<=”写成了“<”。
    这个错误带来的启示就是在写循环条件时一定要特别注意对边界的检查。否则就会失之毫厘谬之千里。
2.        其次，这行代码的语意是每次循环都调用一下sqrt()函数。但实际上，由于sqrt(100)是个常量表达式，所以只调用一次就足够了，多次调用显然是毫无意义没有必要的，因而也是效率低下的。
当然可能有人会辩解说有的编译器在有些条件可以对此进行优化，但这并不表明代码写得不垃圾。自己把代码写的稀里糊涂莫名其妙效率低下，然后把希望寄托给编译器的优化是无论如何都说不过去的。
3.        最后表达式i<sqrt(100)并非绝对地等价于i<10。原因是sqrt()的函数原型是
            double sqrt(double);
    也就是说它的返回值是double类型。由于i作为int类型无法与一个double类型的数据做“<”运算，根据C语言的隐式类型转换规则，i<sqrt(100)这个表达式的真正含义其实是
            (double) i < sqrt(100.)
    亦即是两个double类型进行“<”运算。在C语言中，double等浮点类型并不是用来精确地表达一个值的（虽然个别情况可能也精确），而是用以表示一个实数的近似值的。因而对于上式来说，你即不能指望(double) i准确地表示了i的值，更不能指望sqrt(100.)的值就是准确的10（除非你清楚两件事，第一你的编译器采用的浮点编码制，其次你了解sqrt()这个库函数采用的算法），sqrt(100.)只是10的一个近似值，sqrt(100.)的返回值比10稍微小一点点或大一点点都有可能。因此(double) i < sqrt(100.)这个表示对两个近似值进行近似比较的表达式和i<10这个对两个精确值进行精确比较的表达式的含义完全不同，更不要说这个表达式中的“<”原本应该是“<=”了。
    除此之外

for(i=2;i<sqrt(100);i++)
     for(j=i+1;j<=100;j++)
        {if(a[ i ]!=0&&a[j]!=0)
            if(a[j]%a[ i ]==0)
              a[j]=0;             //把非素数“挖掉”
     }

这段代码，根据前面对筛法的介绍，根本就不是筛法。因为筛法的核心思想是每确认一个素数就“筛”去这个素数的倍数。而“筛”去该素数的倍数，即不需要通过“j++”对数组中该素数后面的数逐个检查，同样不需要进行求余运算（a[j]%a[ i ]），因为只根据下标就可以判断出是否是某个素数的倍数。譬如，若a[2]的值为2是一素数，那么a[2+ a[2]]、a[2+ a[2]+ a[2]]……就一定是2的倍数。所以，即使按照原来不合理的数据结构(int a[101] ; )，这段代码也应该这样写：

for(i=2 ; i*i<=100;i++)
      if(a[ i ]!=0)
         for(j=i+a[ i ];j<=100;j+=a[ i ])
              a[j]=0;

    这才称得上是筛法。可以看到这里 j 并不是每次加1而是加a[ i ]，而且不需要费时的“%”运算，在效率方面与原来的写法天壤之别。这就是筛法的精髓，并不需要做复杂费时的求余运算，只需要做简单的加法就足够了。
样本代码最后一部分是输出素数，这部分的主要缺点是写的不够简洁，从美学的角度看把换行这个很次要的功能写在了一个比较显要的位置也不恰当。这段代码完全可以这样写

for(i=1,n=0;i<=100;i++)
     if( a[ i ] !=0 )                 
      printf( "%5d%s", a[ i ] , ++n%10==0?"\n":""  ) ;  

    下面是正确运用筛法的完整代码：

#define UPPER 100
#define LOWER 1

#define SIZE  ( UPPER - LOWER + 1 )  //数组尺寸 

#define NOT_A_PRIME   0  

#define NUM_PER_LINE  10

#define FOLLOW_BY(n) (( ++ (n) % NUM_PER_LINE == 0 ) ? "\n" : "" )

int main( void )
{
  int num [ SIZE ] , i , n ;
  
  for ( i =0 ; i < SIZE ; i ++ )  //将1~100写入数组 
     num [ i ] = i +1 ; 
  
  num[ 1-1 ] = NOT_A_PRIME ;     //1不是素数,存储1数组元素的下标为1-1 
  
  for ( i =0 ; i * i  <= UPPER ; i ++ )
     if(  num [ i ] !=  NOT_A_PRIME )       //如果num[  i  ]是素数 
     {
        int j ;
        for ( j = i + num [  i  ] ; j < SIZE ; j +=  num [  i  ]  ) 
            num [ j ] = NOT_A_PRIME ;   //num[ i ]的倍数不是素数 
     }      

  for ( i =0 , n =0 ; i < SIZE ; i ++ )
     if( num [  i  ] !=  NOT_A_PRIME )
       printf("%5d%s" , num [  i  ] , FOLLOW_BY(n) ) ;    
      
  putchar('\n');

  return0;
}

    回过头来再来说一下不需要存储1~100这100个数的问题。由于数组下标是连续的整数，实际上程序可以利用这个性质记住各个整数是否是素数就足够了。亦即可以考虑这样的解决方案，用num[0]记录1是否是素数，num[1]记录2是否是素数，……。记录是否是素数并不需要int类型的数据，char类型就足以胜任。这种方案无疑比存储1~100这100个数要更加节省内存。

#include <stdio.h>

#define UPPER 100
#define LOWER 1
#define SIZE  ( UPPER - LOWER + 1 )  //数组尺寸 

#define NOT_A_PRIME   0  
#define IS_A_PRIME    1
  
#define NUM_PER_LINE  10

#define FOLLOW_BY(n) (( ++ (n) % NUM_PER_LINE == 0 ) ? "\n" : "" )

int main( void )
{
  char num [ SIZE ] , i , n ;
  
  for ( i =0 ; i < SIZE ; i ++ )  //先假定1~100都是素数，之后再筛去非素数 
     num [  i  ] = IS_A_PRIME ; 
  
  num[ 1-1 ] = NOT_A_PRIME ;//1不是素数,记录1是否为素数数组元素的下标为1-1 
  
  for ( i =0 ; i * i  <= UPPER ; i ++ )
     if(  num [  i  ] == IS_A_PRIME )       //如果( i + 1 )是素数 
     {
        int j ;
        for ( j = i + ( i +1 ) ; j < SIZE ; j +=  ( i +1 ) ) 
            num [ j ] = NOT_A_PRIME ;   //( i + 1 )的倍数不是素数 
     }      

  for ( i =0 , n =0 ; i < SIZE ; i ++ )
     if( num [  i  ] == IS_A_PRIME )
       printf("%5d%s" , i +1 , FOLLOW_BY(n) ) ;    
      
  putchar('\n');

  return0;
}

    如果希望进一步节约内存，甚至可以考虑用1位来记录某个整数是否为素数。但这已经属于题外话了。
    最后再简略地谈一下求平方根的问题。在前面代码中用i * i  <= UPPER代替了不正确的i<sqrt(100)或i<=sqrt(100)，但这个写法还可以进一步改进。
    由于前n个连续奇数的和恰好是n的平方，所以可以用UPPER依次减去连续的奇数，一旦差小于0，则表明i已经大于√n。