利用辛達拉姆篩進行素數判定

標籤：素數判斷   埃氏篩法   辛達拉姆篩   

1 常規判定方法

       素數判定問題就是對給定的正整數n判定是否為素數。所謂素數，是指恰好有2個約數的整數。因為n的約數都不超過n，所以只需要檢查2~n-1的所有整數是否整除n就能判定是不是素數。不過，我們還能進一步最佳化。如果d是n的約數，那麼n/d也是n的約數。由n=d*n/d可知min(d,n/d)，所以只需要檢查2~的所有整數就足夠了。此時，素數判定的複雜度為O()。代碼實現如下：

int classic(long long a){    long long half=(long long)sqrt(a);    for(int i=2;i<=half;i++)    {        if(a%i==0)  return 0;    }    return 1;}

常規的素數判定思想和實現都非常簡單，但是效率較低。現在存在不少更高效的素數判定方法，比如費馬測試、演算法、數域篩法等，但是這些演算法要麼不精確要麼需要的數學知識過於複雜，在此不做過多介紹。下面我們介紹一種只需要初等數學即可明白的素數判定方法，該方法在任何情況下都會獲得比常規方法好一倍的效能。

2 埃氏篩法

       在介紹新方法之前，我們先回顧一下埃氏篩法，這是一個與輾轉相除法一樣古老的演算法。

       大約在公元前3世紀，古希臘數學家埃拉托色尼提出了一種編造素數表的方法（如表一）。這種方法類似於篩東西，把不要的篩掉，把需要的留下。具體做法是：將從2到N的自然數，按順序排列成2,3,4,5,…,N，然後留下第一個2，划去所有2的倍數；2之後沒划去的第一個數是3，留下3，划去所有3的倍數；在3後面沒劃掉的第一個數是5，留下5，划去所有5的倍數；如此繼續，直至上述一列數中再也沒有可劃的數為止，留下來的便是N以內的一切素數。

表一 埃氏篩法

       如果只對一個整數進行素性判定，常規判定方法或者後面介紹的方法就足夠了。但如果要對許多整數進行素性判定，就需要採用埃氏篩法。以下枚舉演算法摘自《挑戰程式設計競賽》一書：

int prime[MAX_N];bool is_prime[MAX_N];//返回n以內的素數個數int sieve(int n) {int p=0;for(int i=0;i<=n;i++) is_prime[i]=true;is_prime[0]=is_prime[1]=false;for(int i=2;i<=n;i++) {if(is_prime[i]) {prime[p++]=i;for(int j=2*i;j<=n;j+=i) is_prime[j]=false;}}return p;}

3 辛達拉姆篩

       公元1934年，一名年輕的東印度數學生辛達拉姆，提出了一種與埃拉托色尼迥然不同的篩法。它首先列出了一張表，如表二。表的第一行和第一列都是首項為4，公差為3的等差數列。從第二行開始，以後各行也是等差數列，公差分別為5,7,9,11,13……。可以看出，該表其實是一個對稱矩陣。

       辛達拉姆指出： 如果N出現在表中，則2N+1是合數；若N不在表中，則2N+1是素數。證明相當精彩！

表二 辛達拉姆篩構造的表

       先證明前半部分。首先，他寫出來第m行的第一個數：

4+(m-1)×3=3m+1

注意到該行是公差為2m+1的等差數列，所以此行第n列的數是：

3m+1+(n-1)*(2m+1)=2mn+m+n

即第m行n列的數是2mn+m+n。於是2N+1=4mn+2m+2n+1=(2m+1)(2n+1)是合數。

       再證後半部分。要想正面證明2N+1是素數是相當困難的。如果換成等價的逆否命題，即證“若2N+1不是素數，則N必在表中”似乎要容易得多。事實上，若

       2N+1=x*y(x ,y)為整數

       則因2N+1為奇數，x、y也必為奇數。不妨設：

       x=2p+1;y=2q+1

       從而2N+1=(2p+1)(2q+1)=4pq+2p+2q+1。由此可以得到

       N=2pq+p+q。

也即N是表中第p行第q列的數。

       綜上所述，我們證明了辛達拉姆篩法的正確性。例如18不在表中，則2*18+1=37為素數。相反，71在表中，則2*71+1=143是合數。後半部分的證明是反證法的絕佳執行個體。

4 基於辛達拉姆篩的素數判斷   

       有了辛達拉姆篩法，我們該如何判定一個正整數n是否為素數呢？我們只需要判斷(n-1)/2是否在表中即可。由於需要判定的數肯定是奇數（偶數肯定都是合數），所以(n-1)/2肯定可以整除。判斷是否在表中，也即判斷是否存在p和q使得2pq+p+q=(n-1)/2。按照這個思路實現的代碼如下：

int xindalamu(long long a){if(a%2==0) return 0;    long long aa=(a-1)>>1;    long long half=(long long)sqrt(aa/2);    for(int i=1;i<=half;i++)    {        if((aa-i)%(2*i+1)==0) return 0;    }    return 1;}

下面詳細解釋一下上面的代碼。首先判斷使用者輸入是否是偶數，如果是直接返回false。然後將使用者輸入減一除以二賦值給aa，後面就需要判斷aa是否在表中。判斷的方法和常規判斷的方法類似，也是需要驗證某一行是否出現整除的情況，不過這裡公式稍有不同。從上述代碼可以看到，在對aa開方的時候，我們將其除以二。這是為什麼呢？判斷aa是否出現在表中，也即判斷是否存在p和q使得2pq+p+q=aa。這和常規方法中的判斷aa=p*q稍有不同。忽略掉兩個小項，我們需要判斷2pq是否等於aa，也即pq是否等於aa/2，這就是除以2的原因。for迴圈從第一行開始檢查，判斷該行是否存在列滿足2pq+p+q=aa。當行號固定之後，判斷的方法是看aa-p是否整除2p+1（稍作變化即可看出）。

       我們分析一下基於辛達拉姆篩的素數判定複雜度。由for迴圈可以看出，複雜度也是O()，不過常係數比常規方法要小。這是因為，在求aa時有個除以二的操作，在開方時又有一個除以二的操作，所以相當於對a/4進行開方，因而最多隻需要對/2的行進行判斷。而常規方法最多則需要驗證個不同的因子，所以可以看出基於辛達拉姆篩的素數判定方法比常規方法快一倍。

       但是，還有一個疑問需要解決，如果使用者輸入的是一個合數，基於辛達拉姆篩的素數判定方法是否還優於常規判定方法。可以證明，不管輸入的是合數還是素數，基於辛達拉姆篩的素數判定方法始終比常規方法快一倍。證明如下：

       如果使用者輸入的是素數，在分析複雜度的時候已經證明。如果使用者輸入的是合數，則常規方法在遇到a的最小質因子b時返回，而辛達拉姆判定方法會在遍曆到第(b-1)/2行時返回。假設a=b*y，因為a是合數，則肯定存在p行和q列使得2pq+p+q=(a-1)/2。變換得(2p+1)(2q+1)=a。由於辛達拉姆判定方法在遍曆時是從第一行開始的，所以我們可以認為b=2p+1，從而辛達拉姆判定方法在第(b-1)/2行時返回，這也就證明了當輸入是合數（奇數）時，基於辛達拉姆篩的素數判定方法也比常規方法快一倍。

5 基於辛達拉姆篩的素數篩選

       辛達拉姆篩最初的目的和埃拉托色尼篩一樣，也是篩選所有的素數，上面的素數判定算是一個變形應用，篩選所有素數才是篩的最重要目的。不過和埃氏篩法相比，辛達拉姆篩卻不那麼容易獲得所有的素數。從辛達拉姆得到的結論：如果N出現在表中，則2N+1是合數；若N不在表中，則2N+1是素數，我們可以很容易知道不在表中的所有數，其二倍加一包含了所有的奇素數。例如，表中出現的前10個數為4,7,10,12,13,16,17,19,22,24，不在表中的數即為1,2,3,5,6,8,9,11,14,15,18,20,21,23，這些數的二倍加一分別為：3,5,7,11,13,17,19,23,29,31,37,41,43,47，也即前14個奇素數。雖然我們只需要找到所有不在表中的所有數即可獲得所有的奇素數，但是如何找到卻是個麻煩問題。問題在於如何按序遍曆在表中的所有數。本人暫時沒有思路，還請高手幫忙解答，在此謝過。

參考       對於辛達拉姆篩的瞭解是從高中一本數學習題集中獲得的，當時讀了練習後面的拓展閱讀，印象非常深刻，以至於到現在我都留著那本習題集，感謝它。