8. 篩法概述

全名為 埃拉托斯特尼(Eratosthenes，古希臘數學家) 篩法，這裡只概述其大致概念。

假設要找出 1~40(別太大，不然不好說明)內所有質數，先從  1 寫到 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

接著一個數一個數判斷，

(1) 1 非質數也非合數，直接先畫掉
(2) 2 沒被畫掉，是質數，把 2 的倍數全都劃掉

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

(3) 3 沒被劃掉，是質數，把3的倍數也全劃掉

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

(4) 4 被劃掉了，不理它，繼續往下看

(5) 5 沒被劃掉，它是質數，再把 5 的倍數全都劃掉

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

(6) 6已被劃掉，不理它，繼續看

(7) 7 沒被劃掉，它是質數，再把 7 的倍數全劃掉。

(8) 由於 40 的平方根介於 6 與 7 之間，所以上面動作做到 7 即可，接下來剩下的數全都是質數。

看得出來有 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37。

  

9. 篩法 Coding 注意事項

篩法最大優勢是在於找一定數以內的所有質數，在 coding 時要注意，

 動態陣列應用 unsigned char / char / bool，別用 int

上面三個資料型態都吃 1 byte，但 int 卻吃 4 bytes，在做篩法時之 table 只是要紀錄 "是" 或 "不是" 質數，

因此沒必要開那麼大的空間來吃記憶體空間，這樣只會使得 N 之上限值沒那麼大！

( coder 於 OS 可配置之記憶體空間僅有 user space ，XP 為 2GB / 3GB)

這部份單純之 coding 如下所示

void* se_basic(int Max, int* cnt)
{
    // total size: Max+1; un size: 0, 1;Max+1-2 = Max-1
    register int PrimeCnt=Max-1; 
    register int i, j, end;
    byte* p=(byte*)malloc(Max+1);  

    if(p==NULL){
       *cnt=0;
       return NULL;
    }

    memset(p, 1, Max+1),p[0]=p[1]=0;
    end = (int)ceil(sqrt((double)Max));      

    for(i=2; i<=end; ++i){
       if(p[i])
          for(j=i+i; j<=Max; j+=i) PrimeCnt-=p[j], p[j]=0;
    }
    *cnt = PrimeCnt;
    return p;
}

 

上述是最基本型的篩法。

10. 套用 2n /6n 法則

即使是篩法也可應用 2n / 6n 法則。即 2 / 2.3 之倍數都不去判斷，試想一下，2 這個數的倍數就佔了整張表之一半，一開始要刪時便花了近一半的時間，不合效率，於是將 2n 法則套近篩法裡面。

// ==================================
// check odd only
void* se_2n(int Max, int* cnt)
{
    // total size: Max+1; un size: 0, 1; 2 is prime; Max+1-2+1 = Max 
    register int PrimeCnt= (Max>>1)+(Max&1);
    register int i, j, end;
    byte* p=(byte*)malloc(Max+1);  

    if(p==NULL){
       *cnt=0;
       return NULL;
    }

    memset(p, 1, Max+1),p[0]=p[1]=0;
    end = (int)ceil(sqrt((double)Max));

    for(i=3; i<=end; i+=2){
       if(p[i])
          for(j=3*i; j<=Max; j+=(i<<1)) PrimeCnt-=p[j], p[j]=0;
    }
    *cnt = PrimeCnt;
    return p;
}

 

2n 可以套樣，一樣的，6n 、30n 、210n 也可以套用，其實用到 6n 已夠了，210n 實在太麻煩且也沒比較快，30n code 一起附上

void* se_sqrt_6n(int Max, int* cnt)
{
    // total size: Max+1; un size: 0, 1; 2 is prime; Max+1-2+1 = Max 
    register int PrimeCnt;
    register int i, j,  end, gap1;
    byte* p=(byte*)malloc(Max+1);  

    if(p==NULL){
       *cnt=0;
       return NULL;
    }

    memset(p, 1, Max+1),p[0]=p[1]=0;
    end = (int)ceil(sqrt((double)Max));

    for(gap1=4, i=5; i<=end;gap1^=6, i+=gap1){
       if(p[i]){
          for(j=i*i; j<=Max; j+=i+i) p[j]=0;
       }
    }

    for(i=5, gap1=4, PrimeCnt=2; i<=Max; gap1^=6, i+=gap1 ){
       PrimeCnt += p[i];
    }
    *cnt = PrimeCnt;
    return p;
}

void* se_sqrt_30n(int Max, int* cnt)
{
    // total size: Max+1; un size: 0, 1; 2 is prime; Max+1-2+1 = Max 
    register int PrimeCnt;
    register int i, j, k, end;
    byte* p=(byte*)malloc(Max+1);  

    uint a_s[] = {6,4,2,4,2,4,6,2};//step
    uint a_bound = sizeof(a_s)/sizeof(a_s[0])-1;

    if(p==NULL){
       *cnt=0;
       return NULL;
    }

    memset(p, 1, Max+1),p[0]=p[1]=0;
    end = (int)ceil(sqrt((double)Max));

    for(k=0, i=7; i<=end; i+=a_s[k]){
       if(p[i]){
          for(j=i*i; j<=Max; j+=i+i) p[j]=0;
       }
       (k==a_bound)?(k=0): (++k);
    }

    for(PrimeCnt=3, i=7, k=0; i<=Max; i+=a_s[k]){
       PrimeCnt += p[i];
       (k==a_bound)?(k=0): (++k);
    }
    *cnt = PrimeCnt;
    return p;
}

 

 11. 將不必要空間拿掉

以 2n 為例，事實上有一半空間用不到，於是把這一半空間都拿掉

// ==================================
// using half space for 2n rule
void* se_2n_comp(int Max, int* cnt)
{   
    register int pSize = (Max>>1) + (Max&1);
    register int PrimeCnt = pSize;
    register int i, j, step, end;

    byte* p=(byte*)malloc(PrimeCnt);  
    memset(p, 1, PrimeCnt),p[0]=0;

    i = (int)(ceil(sqrt((double)(Max))));
    end = (i - (!(i&1))) >> 1;

    for(i=1; i<=end; ++i){
       if(p[i]){

          // number = (2*i+1)
          // number^2 = (2*i+1) * (2*i+1) = 4*i*i+4*i+1;
          // number^2 of index = number/2 = 2*i*i + 2*i = 2*i*(i+1)
          j = i*(i+1)<<1;
          
          // index of step 
          step = (i<<1)+1;
          for(; j< pSize; j+=step) {
             PrimeCnt-=p[j], p[j]=0;
          }
       }
    }
    *cnt=PrimeCnt;
    // show table
    return p;
}

 

12. 以 bitwise 方式操作

上述已拿掉一半之空間了，但仍不夠！要存一個數是不是質數，只要用 1 bit 去存便可，根本就用不到 1 bytes，這樣只是浪費其它 7 bits 而已，以 bitwise 操作較隱晦，下列程式碼以 2n 為例作參考

#define IndexToNum2(index )(1+((index)<<1))
#define NumToIndex2(number)((number-1)>>1 )
#define SET_BIT(n, bit)((n)|=(0x01 <<(bit)))
#define CLR_BIT(n, bit)((n)&=~(0x01 <<(bit)))
#define CHK_BIT(n, bit)(((n)>>(bit))&1)


void* se_2n_Hcomp(int Max, int* cnt)
{
    // seven register variable.
    register int i, j, Check_end= (Max>>1) + (Max&1); /* check bound to resize */
    register int pSize = (Max >> 6) + 1 - ( (Max&63)==0);
    register int PrimeCnt = Check_end;
    register int end, step;

    uint* p = (uint*)malloc(sizeof(uint)*pSize);
    memset(p, 0xffffffff, sizeof(uint)*pSize), p[0]=0xfffffffe;

    i = (int)ceil(sqrt((double)Max)); /* cal sqrt(N) */
    end = NumToIndex2(i);  /* cal the last check index */ 
    for(i=1; i<=end; ++i){
       if( CHK_BIT(p[i>>5], (i&31))) {
          j = i*(i+1)<<1;
          step = (i<<1)+1;          
          for(; j<Check_end; j+=step){
             PrimeCnt -= ( CHK_BIT(p[j>>5], j&31));
             CLR_BIT(p[j>>5], j&31);
          }
          j-=step;
         // printf("%d\n", IndexToNum2(j)); // check bound.
       }      
    }
    *cnt = PrimeCnt;
    return p;
}

 這些都是經由 bit 運算轉換去推出來的寫法，同時這種操作方式也會使得速度加快，而且快不少。

以 6n 方式進行 bitwise 操作，這部份實不容易，但網路上仍可找到這部份文章說明，於此便不贅述。

 

* 線性篩法

 

上述的篩法是套用了 2n 法則，所有 2 之倍數都不去砍，減少 polling 時間；同樣也有 6n 法則，所有 2 與 3 之倍數都不去砍，又減少了 polling 時間，但這二個方法都會有重覆刪除的現象。線性篩法優點便可保證每個要刪除的數字都只刪除一次，速度為線性，故稱線性篩法。但缺點是需要額外的記憶體空間存質數表。換而言之，如果可定址空間只在 2GB，且欲篩選範圍甚大，此時便不適用線性篩法，因質數表所需空間甚大。 這次以反邏輯方式撰之，基礎型的線性篩法大致如下

#include <stdio.h>
#include <math.h>

#define N 100
int Prime[N]={0};
char p[N]={0};

int main()
{
    int i, j, end, Cnt=0;
    p[0]=p[1]=1;

    for(i=2; i!=N; ++i){
        if(p[i]==0) Prime[Cnt++]=i;
        for(j=0; j!=Cnt && i*Prime[j]<=N; ++j) 
            p[i*Prime[j]]=1;
    }
    
    Cnt=0;
    for(i=0; i!=N; ++i) if(p[i]==0) printf("%d", i), ++Cnt;
    printf("\nCnt: %d\n", Cnt);
    return 0;
}

線性篩法過程中，甚至還可以再套用 2n 、6n 法則

#include <stdio.h>

#define N 100
int Prime[N]={0};
char p[N]={0};

int main()
{
    int i, j, end, gap, Cnt;
    p[0]=p[1]=p[4]=1;

    for(gap=2, i=5; i<N; i+=gap, gap^=6){
        if(p[i]==0) Prime[Cnt++]=i;
        for(j=0; j!=Cnt && i*Prime[j]<=N; ++j) 
            p[i*Prime[j]]=1;
    }
    
    Cnt=0;
    for(i=0; i!=N; ++i) 
        if(i==2 || i==3 || (p[i]==0 && i%2!=0 && i%3!=0)) 
            printf("%d", i), ++Cnt;
    printf("\nCnt: %d\n", Cnt);
    return 0;
}

甚至可再進一步用 bit 存資訊、以 bitwise 操作，於此便不再示範。

參考資料

wiki - Mersenne prime
wiki - Fermat number
wiki - Carmichael number

wiki - Sieve of Eratosthenes

wiki - Modular exponentialiation
wiki - Montgomery reduction

wiki - Fermat primality test
wiki - Miller-Rabin primality test
wiki - Lucas Lehmer test for Mersenne numbers
wiki - AKS primality test

wiki - Prime-counting function
wiki - Euler's totient function

arrow
arrow
    文章標籤
    質數
    全站熱搜
    創作者介紹
    創作者 edisonx 的頭像
    edisonx

    Edison.X. Blog

    edisonx 發表在 痞客邦 留言(0) 人氣()

    E-mail轉寄