2009年10月27日 星期二

質數

 
忙裡偷閒,純粹寫好玩的 :-)

列出數字 1~n 間的質數:
int bPrime (int a, int b) 
{
    return b<2 && printf (" %d", a) || 
           a%b && bPrime (a, b-1);
};
void print (int n) 
{ 
    if (n>2) print (n-1); 
    bPrime (n, n-1); 
}
int main (/*daviddr 091028*/)
{
    print (20);
    return getch();
}

以昇鼏方式判斷是否為質數:
bool bPrime (int n, int i=2) 
{
    return i==n && printf (" %d", n) || 
           n%i  && bPrime (n, i+1);
};
void print (int n) 
{ 
    for (int i=2; i<=n; i++)
        bPrime (i); 
}
int main (/*daviddr 091028*/)
{
    print (20);
    return getch();
}

略去輸出,便為:
int bPrime (int n, int i=2) 
{
    return i==n || n%i && bPrime (n, i+1);
};

「似乎」較快,但潛藏危機的寫法:
(false 必為 0,但 true 不一定是 1)
int bPrime (int n, int i=2) 
{
    return i==n || n%i && bPrime (n, i+(i>2)+1);
};

或者也可以用機械碼來寫,
不過,會被平台綁定 (此例為IA32):
void print_prime (int a, int b) 
{ 
    long long _[] = {
        0x0000ccec81ec8b55, 0xff34bd8d57565300,
        0xb800000033b9ffff, 0x45c7abf3cccccccc,
        0x8b09eb00000002f8, 0xf8458901c083f845,
        0x117d08453bf8458b, 0x85f87df79908458b,
        0xeb04ebc0320475d2, 0xe58b5b5e5f01b0de,
        0xc35d
    };
    for (;a <= b; ++a)
        if (((bool(*)(int))&_)(a)) 
            printf (" %d", a); 
}
int main (/*daviddr 091030*/)
{    
    int a, b;
    printf ("\n輸入數值範圍: ");
    scanf  ("%d %d", &a, &b);
    print_prime (a, b); 
}

當n為1至100間時,可用位元篩來印:
int main ()
{
    int  j, n;
    long long i = 150964650272191;

    printf ("Input: ");
    scanf ("%d", &n);               //輸入數字

    if (2 <= n) putchar ('2'); 
    n = (n-3)/2;            
    
    for (j=0; j<n; j++, i>>=1)
        if (i&1) printf (" %d", j*2+3);
}

n值更大時,可逐漸編織篩子:
int  i, j, p;
bool sieve [5000];                  //儲存 3,5,7,...等奇數, 通式為 2*i+3

int main ()
{
    int n;
    scanf ("%d", &n);               //輸入數字
    if (n > 2) printf("%8d",2);
    n = (n-3)/2;                    //算出 n 對應的 sieve 奇數索引
    for (i=0; i<=n; i++)
        if (!sieve[i]) {            //若此奇數仍未被刪掉,表示他是質數
            p = i+i+3;
            for (j=p+i; j<=n; j+=p)
                sieve[j] = true;    //把 2*i+3 的倍數都刪掉
            printf ("%8d", p);      //印出此質數
        }
}

當篩子材料不夠用時,
可接替代換成基本定義解:
int main()
{
    int i, j, n;
    scanf ("%d", &n);                 //輸入數字
    for (i=2; i<=n; ++i) {
        for (j=2; j<i; ++j)     
            if (i%j == 0) goto Next;  //若不是質數就換下一個
        printf("%8d",i);
    Next:;
    }
}

因一數的最大因數,必 <= 其根號值。
於是可用如下方式跳過不需計算的數:
int main (/* [式1] */)
{
    int i, j, n, sn;

    printf ("輸入大於 3 的數: ");
    scanf  ("%d", &n);
    printf ("%8d%8d", 2, 3);

    for (i=4; i<=n; i++) {
        if (i%2==0 || i%3==0) continue;
        sn = sqrt ((double)i);
        for (j=5; j<=sn; j+=6)
            if (i%j==0 || i%(j+2)==0) 
                goto Next;        
        printf ("%8d", i);
    Next:;           
    }
}

此法可以不斷擴展下去,譬如:
int main (/* [式2] */)
{
    int i, j, n, sn;

    printf ("輸入大於 7 的數: ");
    scanf  ("%d", &n);
    printf ("%8d%8d%8d%8d", 2, 3, 5, 7);

    for (i=8; i<=n; i++) {
        if (!(i%2 && i%3 && i%5 && i%7)) 
            continue;            
        sn = sqrt ((double)i);
        for (j=7; j<=sn; j+=30)
            if (!(i%j      && i%(j+ 4) &&  
                  i%(j+6)  && i%(j+10) &&
                  i%(j+12) && i%(j+16) &&
                  i%(j+22) && i%(j+24)))
                goto Next;        
        printf ("%8d", i);
    Next:;           
    }   
}

[式2] 將 [式1] 的計算量由 1/3 減至 4/15,
再推展下去,減低的計算量會趨近一個底線...
恩.. 所以推展太多項也不見得好。

若題目不是找1~N裡的所有質數,而是「驗證」質數,
且質數範圍限定為 32-bits int,如:-2147483648 ~ 2147483647
或正整數 0 ~ 4294967295,則查表法更顯威力~
因為可以用 Fermat 小定理 + Carmichae number 來測試,
32-bits int range 中,僅有 1118 個 Carmichae number,
先全算出來放入陣列,然後用 binary search 縮小測試數目,
便可達到快速驗證。

若數值範圍為 16-bits,則可直接用 bit 紀錄 sieve 篩,
達到小量 table 儲存空間、O(1) 存取速度的質數檢測。

若是 32-bits 以上的大數,用查表法則不一定高效,
因為當 memory 配置過大,cache 失誤率也將提高...
對於這種大數,可用 Rabbin-Miller 做隨機(確率)驗證,
如 Java 的 BigInteger.isProbablePrime ;
或用 Agrawal-Kayal-Saxena 做精確驗證,
但 AKS 速度慢,雖然後來 Bernstein 有為它做加速...

沒有留言:

張貼留言