質數的定義
一個數，如果只有1和它本身兩個因數，這樣的數叫做質數，又稱素數。
在上文 《素數算法大全，及C程序實現優化詳解 (一) 試除法》中我們已經探討了求解素數的一類算法，并且將試除法從最初的低效版本優化的高效的V2。那么，還有沒有其它更佳算法呢？這就是下面三藏要和大家探討的內容
合數過濾篩選法
算法描述：我們知道，素數N不能被2~(N-1)間的任何數整除；反過來看，只要能被2~(N-1)間的任何數整除的N，都不是素數。所以我們可以采用一個簡單的排除法：就是對N以內的所有數，只要逐個去除值為2~(N-1)的倍數的數，剩下的就是素數。
C語言實現
// 合數過濾篩選法 Ver1
// 參數：n 求解n以內(包括n)的素數
// 返回值：n以內素數個數
int CompositeNumFilterV1(int n)
{
int i, j;
// 素數數量統計
int count = 0;
// 分配素數標記空間，結合后文思考為何+1
char* flag = (char*)malloc( n+1 );
// 初始化素數標記
for (i=2; i<=n; i++)
{
// 為什么*(p+i)要寫成flag[i]呢？可讀性更佳爾
flag[i] = 1;
}
// 寫程序要注意排版和留空，方便閱讀，也可減少出錯幾率
// 以2~(N-1)為因子過濾合數
for (i=2; i < n; i++)
{
for (j=2; i*j <= n; j++)
{
// i*j是由i,j兩整數相乘而得，顯然不是素數
flag[i*j] = 0;
}
}
// 統計素數個數
for (i=2; i<=n; i++)
{
// 其實if(flag)就其同樣作用了，但這么寫是有留言的
// 請參閱《C語言程序設計常見錯誤剖析及解決之道》一文
if (1 == flag[i]) count++;
}
// 因輸出費時，且和算法核心相關不大，故略
// 釋放內存，別忘了傳說中的內存泄漏
free(flag);
return count;
}
在上文給出的main函數中以不同參數調用CompositeNumFilterV1函數，得到執行結果如下：
[100000]以內素數個數：9592, 計算用時：15毫秒
[1000000]以內素數個數：78498, 計算用時：125毫秒
[5000000]以內素數個數：348513, 計算用時：2578毫秒
[10000000]以內素數個數：664579, 計算用時：6281毫秒
注：因程序是非獨占性運行的，所以時間不是完全精確的，但基本能反映實情
顯然，比上文中的試除法要快，而且誰都可以看到上例是一個未經優化的粗陋版本，好多地方是三藏故意采用比較低效做法，為了與后文的優化版比較，凸顯優化之重要，也為了初學者記住別采用類似低效做法，下面我們開始優化之旅
優化分析
上面CompositeNumFilterV1函數存在的問題有：
- 在外層循環，需要一直執行到n-1嗎？不要，因為n/2~n-1間的數顯然不能整出n
- 在內層循環中重復使用i*j顯然是低效的，考慮到計算機中加減運算速度比乘除快，可以考慮變乘法為加法
- 在循環修改flag過程中，其實有很多數會被重復計算若干次，比如6=2*3=3*2，會被重復置0，類似操作很多，所以我們得設法避免或減少flag重復置0
據上述分析，我們可將程序優化如下：
// 合數過濾篩選法 Ver2
// 參數：n 求解n以內(包括n)的素數
// 返回值：n以內素數個數
int CompositeNumFilterV2(int n)
{
int i, j;
// 素數數量統計
int count = 0;
// 分配素數標記空間，明白+1原因了吧，因為浪費了一個flag[0]
char* flag = (char*)malloc( n+1 );
// 初始化素數標記，要高效點咯
flag[2] = 1;
// 注意是i<n不是上例中的i<=n了，理由自思
for (i=3; i<n; i++)
{
flag[i++] = 1;
// 偶數自然不是素數，直接置0好了
flag[i] = 0;
}
// n為奇數
if (n%2 != 0)
{
flag[n] = 1;
}
// 從3開始filter，因為2的倍數早在初始化時代就干掉了
// 到n/2止的理由還要說嗎
for (i=3; i <= n/2; i++)
{
// i是合數，請歇著吧，因為您的工作早有您的質因子代勞了
if (0 == flag[i]) continue;
// 從i的2倍開始過濾，變乘法為加法
for (j=i+i; j <= n; j+=i)
{
flag[j] = 0;
}
}
// 統計素數個數
for (i=2; i<=n; i++)
{
if (flag[i]) count++;
}
// 因輸出費時，且和算法核心相關不大，故略
// 釋放內存，別忘了傳說中的內存泄漏
free(flag);
return count;
}
再來調用CompositeNumFilterV2得到執行結果：
[100000]以內素數個數：9592, 計算用時：n太小，時間精度不夠
[1000000]以內素數個數：78498, 計算用時：31毫秒
[5000000]以內素數個數：348513, 計算用時：453毫秒
[10000000]以內素數個數：664579, 計算用時：1062毫秒
[100000000]以內素數個數：5761455, 計算用時：12973毫秒
哇哇，比昨天的試除發快了好多倍，可見算法的威力，值得好好學習，別說學算法沒用咯。
上例著那個計算一億以內的素數只要約13秒，應該算不錯了，今天是否可以休息了呢？No，我們要追求極限！
int CompositeNumFilterV3(int n)
{
int i, j;
// 素數數量統計
int count = 0;
// 分配素數標記空間，明白+1原因了吧，因為浪費了一個flag[0]
char* flag = (char*)malloc( n+1 );
// 干嘛用的，請仔細研究下文
int mpLen = 2*3*5*7*11*13;
char magicPattern[mpLen];
// 奇怪的代碼，why，思考無法代勞，想！
for (i=0; i<mpLen; i++)
{
magicPattern[i++] = 1;
magicPattern[i++] = 0;
magicPattern[i++] = 0;
magicPattern[i++] = 0;
magicPattern[i++] = 1;
magicPattern[i] = 0;
}
for (i=4; i<=mpLen; i+=5)
magicPattern[i] = 0;
for (i=6; i<=mpLen; i+=7)
magicPattern[i] = 0;
for (i=10; i<=mpLen; i+=11)
magicPattern[i] = 0;
for (i=12; i<=mpLen; i+=13)
magicPattern[i] = 0;
// 新的初始化方法,將2,3,5,7,11,13的倍數全干掉
// 而且采用memcpy以mpLen長的magicPattern來批量處理
int remainder = n%mpLen;
char* p = flag+1;
char* pstop = p+n-remainder;
while (p < pstop)
{
memcpy(p, magicPattern, mpLen);
p += mpLen;
}
if (remainder > 0)
{
memcpy(p, magicPattern, remainder);
}
flag[2] = 1;
flag[3] = 1;
flag[5] = 1;
flag[7] = 1;
flag[11] = 1;
flag[13] = 1;
// 從17開始filter，因為2,3,5,7,11,13的倍數早被kill了
// 到n/13止的，哈哈，少了好多吧
int stop = n/13;
for (i=17; i <= stop; i++)
{
// i是合數，請歇著吧，因為您的工作早有您的質因子代勞了
if (0 == flag[i]) continue;
// 從i的17倍開始過濾
int step = i*2;
for (j=i*17; j <= n; j+=step)
{
flag[j] = 0;
}
}
// 統計素數個數
for (i=2; i<=n; i++)
{
if (flag[i]) count++;
}
// 因輸出費時，且和算法核心相關不大，故略
// 釋放內存，別忘了傳說中的內存泄漏
free(flag);
return count;
}
再看CompositeNumFilterV3執行結果：
[1000000]以內素數個數：78498, 計算用時：15毫秒
[5000000]以內素數個數：348513, 計算用時：203毫秒
[10000000]以內素數個數：664579, 計算用時：515毫秒
[100000000]以內素數個數：5761455, 計算用時：6421毫秒
再次優化后速度提升了又一倍左右，三藏不禁有點滿足了，睡覺也！