[C Programming] 線性篩法

在做篩法求質數的問題時，我們知道在刪除非質數的資料時，
有很多是重複刪除的；比如說如果有一個數是3*7*17*23，
那麼在刪除3的倍數時會刪到它，刪除7、17、23的倍數時也會刪它，
當質數範圍很大的時候，重複的狀況就會越來越多，
因而程式效率就會大打折扣。
基於篩法的觀念，請寫一個程式，在刪除非質數時「絕對」不做重複的工作。

解答：
我們使用所謂的線性篩法，
基本想法是先找第一個質數(2)，刪除2^2、2^3、2^4...
接著刪除2*3、2^2 *3、2^3 * 3、.......
接著是2*5、2^2 * 5、2^3 * 5、......以此類推
每次我們都先找一個質數p，接著找下一個比這個質數大，
但還沒有被刪去的數q(注意q不一定是質數!)，
刪去p^i *q (i=1,2,...)。由於每次刪去的數都為當時的q*(質數p的i次方)，
所以當找到下一個沒被刪去的q'時，不存在q* p^i = q' * p^i'的可能性。
因為假如等號成立，q' = q* p^(i-i')，i != i'，所以q' = q*(p^j)，
如此一來，在以q作刪去的時候就會刪掉q'了，矛盾。
依此推論，這個篩法不會去篩重複的數。

對於所有合成數來說應該都可以寫成如下的形式：
r = p^i * q，p為r的最小的質因數。
由於我們的q是使用所有大於p且沒有被刪去的數，
也就是說對於所有含有質因數p在內的數，應該都會被篩掉。
因為從頭開始篩到x，因為不論如何這之前的數都會進篩選中，
也就是所有x之前的質因數都被選為p給篩過去了，
所以小於x的質數全部都進行篩選過，
那麼現在x只能是質數，否則它就會被篩掉。
換句話說，這樣一個一個掃過去，選到的p必定是質數。

這是一個很不嚴謹的說明而非證明，
只是嘗試將線性篩法的概念表達出來。
想要了解更完備的證明，請參照這篇paper：
A Linear Sieve Algorithm for Finding Prime Numbers
在實作上的部分來說，
要找到2到n的所有質數，那麼p^2最大不能超過n，
同理，p^i*q也不能超過n，
篩選部分的函式大概可以寫成如下的三層迴圈：
for(p=2; p*p<=n ; NEXT(p))
     for(q=p; p*q<=n; NEXT(q))
         for(m=p*q; m<=n; m=p*m)
             REMOVE(m);

我們可以用Double Linked List來做數列的連結和刪去。
這當中往前找上一個值的prev[i]的值為i-1，
往後找下一個值的next[i]的值則為i+1，
prev[2]和next[n]則到頭尾了，設為NULL。
刪除某個數m呢?
我們就把prev[m]的next指向m的下一個，
即next[prev[m]] = next[m]；
next[m]的prev指向m的上一個，
即prev[next[m]] = prev[m]。
最後，要抓下一個待篩的元素，
就令 x = next[x]，就可以得到x的下一個元素了。

Code：

#include <stdio.h>
#include <stdlib.h>

#define MAXSIZE 1000
#define NEXT(x) x = next[x]
#define REMOVE(x) { next[prev[x]]=next[x]; prev[next[x]]=prev[x];}
#define INITIAL(n) { unsigned long i;    \
    for(i=2; i<=n; i++)      \
     prev[i]=i-1, next[i]=i+1; \
    prev[2]=next[n]=NULL;   \
  } 

void L_Sieve(unsigned long n)
{
 unsigned long prev[MAXSIZE+1], next[MAXSIZE+1];
 unsigned long prime, q, i, mult, cnt = 0;
 INITIAL(n);
 for(prime=2; prime*prime <= n; NEXT(prime))
  for(q=prime; prime*q <= n; NEXT(q))
   for(mult=prime*q; mult <= n; mult *= prime)
    REMOVE(mult);
 for(i = 2; i!= NULL; NEXT(i))
 {
  if(cnt % 8 == 0) printf("\n");
  printf("%6ld", i);
  cnt++;
 }
 printf("\n\nTotal Prime Count: %ld \n", cnt);
}

int main(int argc, char *argv[])
{
 L_Sieve(1000);
 system("PAUSE");
 return 0;
}