素數篩法詳解

我們直接抓住素篩的核心：

定理: 如果n不是素數, 則n有滿足1<d<=sqrt(n)的一個"素數"因子d.(若沒有素因子，則n是素數)

證明:

I1. 如果n不是素數, 則n有滿足1<d<=sqrt(n)的一個因子d.

I2. 如果d是素數, 則定理得證, 算法終止.

I3. 令n=d, 并轉到步驟I1.

由于不可能無限分解n的因子, 是以上述證明的算法最終會停止.

代碼:

// primes[i]是遞增的素數序列: 2, 3, 5, 7, ...
// 更準确地說primes[i]序列包含1->sqrt(n)範圍内的所有素數
// 如何高效構造primes[]？需要多少個素數才能判斷2^31-1内的所有素數？一會會有說明


int isPrime(int primes[], int n)
{
    if(n < 2) return 0;


    for(int i = 0; primes[i]*primes[i] <= n; ++i)
        if(n % primes[i] == 0) return 0;//對應上述“n有素因子，n是合數”


    return 1;
}      

構造素數序列primes[i]: 2, 3, 5, 7, ...

由4的算法我們知道, 在素數序列已經被構造的情況下, 判斷n是否為素數效率很高;

但是, 在構造素數序列本身的時候, 是否也可是達到最好的效率呢?

事實上這是可以的! -- 我們在構造的時候《《完全可以利用已經被構造的素數序列!》》

假設我們已經我素數序列: p1, p2, .. pn

現在要判斷pn+1是否是素數, 則需要(1, sqrt(pn+1)]範圍内的所有素數序列,

而這個素數序列顯然已經作為p1, p2, .. pn的一個子集被包含了!

// 構造素數序列primes[]

void makePrimes(int primes[], int num)
{
    int i, j, cnt;
    
    primes[0] = 2;
    primes[1] = 3;
    
    for(i = 5, cnt = 2; cnt < num; i += 2)//cnt為primes下标，i+=2篩去了偶數
    {
        int flag = 1;
        for(j = 0; primes[j]*primes[j] <= i; ++j)
        {
            if(i%primes[j] == 0)
            {
                flag = 0; break;
            }
        }
        if(flag) primes[cnt++] = i;
    }
}      

目前的主流PC中, 一個整數的大小為2^32. 如果需要判斷2^32大小的數是否為素數,

則可能需要測試[2, 2^16]範圍内的所有素數(2^16 == sqrt(2^32)).

在對[2, 2^16]範圍内的素數進行統計, 發現隻有6542個素數:

p_6542: 65521, 65521^2 = 4293001441 < 2^32, (2^32 = 4294967296)

p_6543: 65537, 65537^2 = 4295098369 > 2^32, (2^32 = 4294967296)（在實際運算時unsigned long x = 4295098369;将發生溢出, 為131073）

分析到這裡我們可以看到, 我們隻需要緩沖6543個素數, 我們就可以采用上述算法高效率地判斷[2, 2^32]如此龐大範圍内的素數!