close

這是系列文，前四篇如下

 

[C&++] 深入質數 (1/n) - 試除法單一測試

 
 
 
 
這篇文章主要是裡面議題一些 note ，這些  note 主要是網友的問題，自己又回來摸一下，不過時間有限，沒做太多實作，只把問題點與解決方案稍加紀錄。
 
 
[1] 質數表
 
 
有人問，建立質數表到底有什麼用？
 
答案是它本身是一種可以加快判斷某數 N 是否是質數。為何？試想一下在第一篇文章，最終之「單一試除」，用的是 6n +- 1 去試除，只要試除到 sqrt(N) 即可。可能有個特性沒講清楚。
 
一個合數，一定可以拆成質數之乘積，所以它可以拿來加速判斷一數 N 是否為質數。
 
換句話說，如果我要判斷 91 是不是質數時，步驟如下
 
(1) sqrt(91) = 9.xxxx
(2) 小於等於 9.xxxx  的質數有 { 2,3,5,7 }
(3) 拿上面的那些數去試除
 
 
如果可以除盡時就代表它不是質數了。那質數表到底要建多大？答案是真的不知道、不一定、不確定。但如果真的有需求的話，可以直接寫到檔案裡面去，一直開著做準備。
 
 
另一用途，它可以拿來做較快速之質因數分解。
 
拿剛剛的 273  為例子，拿 {2,3,5,7} 去試除，馬上就分解出來。
 
 
甚至，若一 process 會進行兩次篩法，拿這份質數表 「應」可加快第二次篩法速度，但這應該很少見了，不多提。
 
 
 
 
[2] 怎麼線性篩法好像怪怪的？就算正確也沒比傳統篩法快啊
 
 
 
傳統篩法有個缺點，每個合數可能不只刪除一次，而且可能有些合數會被刪很多次；但線性篩法就保證一定只有刪除一次。但問題又是，我怎麼知道質數表預設到底要設多大？
 
 C++ 用 vector 去存就沒這個問題，C 的話可能比較麻煩，但其實是有一個上限公式存在，我找到的 N >5xx 時，質數上限如下
 
PSIZE = (size_t)(N/lnN * (1.0 + 1.2762 / lnN));
 
這個公式當時我有比較過，比其他公式上限來得小些，比實際個數來得大些。另一份線性篩法 implement 如下
 
Code Snippet
  1. typedef unsigned char byte;
  2.   typedef unsigned intuint;
  3. typedef unsigned long long uint64
  4. byte * SieveLinear(const size_t N, size_t *c)
  5. {
  6.     size_t PSIZE;
  7.     size_t i, j, idx, pn=0, cnt=N-1;
  8.     uint * Prime;
  9.     byte * p;
  10.     double lnx = log((double)N);
  11.  
  12.     PSIZE = (size_t)(N/lnx * (1.0 + 1.2762 / lnx));
  13.  
  14.     Prime = (uint*)malloc(sizeof(uint)*PSIZE);
  15.     if(!Prime) {*c=0; return NULL;}
  16.     p = (byte*)malloc(N+1);
  17.     if(!p){free((void*)Prime); *c=0; return NULL;}
  18.  
  19.     memset(p, 1, N+1), p[0]=p[1]=0;
  20.  
  21.     for(i=2; i<=N; ++i){
  22.         if(p[i]) Prime[pn++] = i;
  23.         for(j=0; j<pn  && (idx = i*Prime[j])<=N ; ++j) {
  24.             cnt-=p[idx], p[idx]=0;
  25.             if(!(i%Prime[j])) break;
  26.         }
  27.     }
  28.     *c = cnt;
  29.     free(Prime);
  30.     return p;
  31. }
  32.  
  33. int SieveLinearCheck(const size_t n, const byte *p)
  34. {
  35.     return p[n];
  36. }
 
上面的 Prime 是質數表，最後被我釋放掉了，程式若還有要用到的話可以先放著別砍。然後我自己在 [C&++] 深入質數 (2/n) - 埃拉托斯特尼篩法 這篇有提到線性篩法可考慮用 2n / 6n 法則，是可以套用沒錯， code 寫出來大致如此。
 
Code Snippet
  1. byte * SieveLinear6n(const size_t N, size_t *c)
  2. {
  3.     size_t PSIZE;
  4.     size_t i, j, idx, gap,pn=0;
  5.     size_t cnt=2+(N+1)/6*2+(N%6==5 || N%6==0);
  6.     uint * Prime;
  7.     byte * p;
  8.     double lnx = log((double)N);
  9.  
  10.     PSIZE = (size_t)(N/lnx * (1.0 + 1.2762 / lnx));
  11.  
  12.     Prime = (uint*)malloc(sizeof(uint)*PSIZE);
  13.     if(!Prime) {*c=0; return NULL;}
  14.     p = (byte*)malloc(N+1);
  15.     if(!p){free((void*)Prime); *c=0; return NULL;}
  16.  
  17.     memset(p, 1, N+1), p[0]=p[1]=0;
  18.  
  19.     for(gap=4,i=5; i<=N; i+=(gap^=6U)){
  20.         if(p[i]) Prime[pn++] = i;
  21.         for(j=0; j<pn  && (idx = i*Prime[j])<=N ; ++j) {
  22.             p[idx]=0;
  23.             cnt-=( (idx%3) && (idx&1U) );
  24.             if(!(i%Prime[j])) break;
  25.         }
  26.     }
  27.     *c = cnt;
  28.     free(Prime);
  29.     return p;
  30. }
 
 
但有個重點蠻多人都忽略的，這段 code 有個地方會有問題，拿去跑 10 億問題就會出來，
 
(idx = i*Prime[j])

這個 idx 可能會導致溢位，這也是常見到的錯誤。原因應不用多說了。解決方式有兩種，第一次是將 idx 宣告成 uint64，並將 i*Prime[j] ，先將 i 轉型成 uint64。但轉成 uint64 之後速度整個拖慢，於是用一點小技巧去改善。

Code Snippet
  1. for(j=0; j<pn  && (idx = i*Prime[j])<=N ; ++j) {
  2.     if(idx < i || idx < Prime[j]) break; // overflow
  3.     p[idx]=0;
  4.     cnt-=( (idx%3) && (idx&1U) );
  5.     if(!(i%Prime[j])) break;
  6. }
 
 
 
 話說回來，就算是線性篩法，也沒普通篩法來得快吧？
 
原因是普通篩法可以用 bitwise 寫法，增加 cache-hit，但線性篩法會不停取用 array element，甚至做二次索引，速度要比普通篩法來得快會有難度 ( 如果普通篩法有寫好的話確實會比線性篩法來得快沒錯)。
 
 
[3] 最快、最神奇的 Code ?
 
 
 
說實話，這段 code 我到現在都還沒能看得懂，也由於我沒能看得懂就不再抄過來了，附上原著的連結 - [原创]求质数 之 筛法(C语言描述) 。
 
該文作者 blog 解釋過，唯程式碼和解釋我實在沒辦法湊在一塊兒。
 
筆者自己寫的是 bitwise-2n 法則，自己環境下10 億跑 10 秒；這份 code 跑 6.5 秒。
 
 
 
另要聲明，這段 code 並不是最快的 code，但無庸致疑想法很好，不過 code 寫得就真的有點... 隱晦了。
 
目前筆者看到最快的 code 有用到 mpi : parallel Sieve Eratosthenes [ updated ] ，筆者的電腦跑 10 億是 3 秒，當然我也不敢保證，這份是目前流傳最快的 code。
 
 
 
[4] 費馬測試 / 拉賓測試是不是有問題啊？？？
 
 
 
我曾介紹過這兩種質數測試方式，[C&++] 深入質數 (3/n) - 質數測試，裡面有含 code。
 
這裡要講，那份 code 只是 present 其概念，如果拿來測 64 位元之質數，通常會出包，
 
這部份網路蠻多也是寫概念，實際跑會有問題，原因又是溢位問題，關鍵在於裡面用到了很多 c = a*b % m 這東西，但 a*b 本身就有溢位之可能性。若改成 c = (a%m) * (b%m) % m，還是可能會在 (a%m) * (b%m) 時發生溢位 ( a%m 及 b%m 都可能大於 32 bits)，
 
故這裡最好自己再寫一個副函式出來做。筆者拿費馬測試做修改
 
Code Snippet
  1.  
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include <time.h>
  5. #include <math.h>
  6.  
  7. typedef unsigned long long uint64;
  8.  
  9. uint64 random(uint64 low, uint64 up)
  10. {
  11.     uint64 x =
  12.         ( (uint64)rand() & 0xFFF0000000000000ULL >> 0 ) |
  13.         ( (uint64)rand() & 0xFFF1000000000000ULL >> 12 ) |
  14.         ( (uint64)rand() & 0xFFF1000000000000ULL >> 25 ) |
  15.         ( (uint64)rand() & 0xFFF1000000000000ULL >> 38 ) |
  16.         ( (uint64)rand() & 0xFFF1000000000000ULL >> 41 ) ;
  17.     return x % (up-low+1) + low;
  18. }
  19.  
  20. // ----------------------------------------------------------
  21. // MultMod(a,b,m) a*b %m 
  22. uint64 MultMod(uint64 a, uint64 b, uint64 m)
  23. {
  24.     uint64 ret=0ULL;
  25.     for(a%=m; b ; b>>=1ULL){
  26.         if(b&1ULL && ((ret+=a)>=m)) ret-=m;
  27.         if((a<<=1)>=m) a-=m;
  28.     }
  29.     return ret;
  30. }
  31.  
  32. // ----------------------------------------------------------
  33. // Montgomery  (a^b % m)
  34. uint64 Montgomery(uint64 a, uint64 b, uint64 m)
  35. {
  36.     uint64 k=1;
  37.     a%=m;
  38.     for(a%=m; b!=1ULL; b>>=1){
  39.         if((b&1)!=0) k= MultMod(a, k, m);// (a*k) % m;
  40.         a = MultMod(a,a,m); //a = (a*a % m);
  41.     }
  42.     return MultMod(a, k, m); // (a*k) % m;
  43. }
  44. // ----------------------------------------------------------
  45. //
  46. int IsPrimeFematTest(uint64 n,uint64 Times)
  47. {
  48.     unsigned i, rnd;
  49.  
  50.     if(n==2 || n==3) return 1;
  51.     else if(n%2==0 || n%3==0 || n%5==0 || n<2) return 0;
  52.  
  53.     if(Times<10U) Times=10U; // Lower Bound
  54.     for(i=0 ; i<Times; ++i){
  55.         rnd = random(2U, n-2);
  56.         if(Montgomery(rnd, n-1, n)!=1ULL)
  57.             return 0;
  58.     }
  59.     return 1;
  60. }
  61.  
  62. int main()
  63. {
  64.     const uint64 N = 6312646216567629137ULL;
  65.     const uint64 T = 100ULL;
  66.     if(IsPrimeFematTest(N, T)) printf("%llu is a primer\n", N);
  67.     else printf("%llu is not a primer\n", N);
  68.  
  69.     getchar();
  70.     return 0;
  71. }
 
 
重點在於，我寫了一份 MultMod 副函式，去做 a*b % m 之動作。它相當低效，若要效率高一些，那就是要實做 128-bits 之大數，這裡不贅述。另米勒拉賓測試改法也類似。最後這部份測試有興趣可再去 K AKS Primer Test。
 
 
[5] 我想確定 64-bits 某個數是否為質數，不要用測試法可以嗎？
 
 
 
嗯，這個蠻特別的，不用測試法實際上還蠻花時間的，只是剛好 64-bits 還算是 "有機會" 可以判斷。
 
假設是 N = 2^50-1，用試除法大概會跑 30 秒。
 
(1)  sqrt(N) = 2^25 ，於是建立 2^25 以內的質數表 (可以用普通篩法沒關係，總之要可以 polling 2^25 以下所有質數)，不到 1 秒。
 
(2) 所以，拿 N 去試除 2^25 以內的所有質數，如果可以整除的話就不是質數了。
 
速度應比試除法快些，但若真的是 2^64-1 時，必須建立 2^32 之質數表，這個篩法就不一定做得出來了。
 
 
 
[6] 我看過一個題目，要求 1000000000000 ~ 1000100000000 之間有幾個質數
 
 
 
注意  1000000000000 ~ 1000100000000，實際區間可調整成 [0:100000000] ，概念上還是用篩法，不過要篩兩次。
 
首先一樣求 sqrt (1000100000000) 以內所有質數，假設叫 Prime[]，然後 拿 Prime 去刪剛調整好的  [0:100000000]，概念不難，但細節要想一下才可，就不再實作了。
 
然後，這種題目必須要先求 sqrt (1000100000000)，這個數值最好事先算好，不然 C 沒辦法精準表達，因 sqrt 裡面引數只有 float / double / long double，轉過去時會造成精度遺失，最後算出來之 sqrt 會偏小。
 
堅持要自己寫 uSqrt 的話下面參考
 
Code Snippet
  1. // ----------------------------------------------------------
  2. uint64 Sqrt64NewtonFloor(uint64 value)
  3. {
  4.     uint64 tmp, v=value>>1;
  5.     if(!v) return value; // 0 return 0, 1 return 1
  6.  
  7.     do{
  8.         tmp=v;
  9.         v = (tmp + value / tmp) >> 1;
  10.     }while( v!=tmp && v!=tmp+1ULL && tmp!=v+1ULL );
  11.     return v - (v*v > value);
  12. }
  13.  
  14. uint64 Sqrt64NewtonCeil(uint64 value)
  15. {
  16.     uint64 tmp, v=value>>1;
  17.     if(!v) return value; // 0 return 0, 1 return 1
  18.     do{
  19.         tmp=v;
  20.         v = (tmp + value / tmp) >> 1;
  21.     }while( v!=tmp && v!=tmp+1ULL && tmp!=v+1ULL );
  22.     return v + (v*v < value);
  23. }
  24. // ----------------------------------------------------------
  25. uint64 Sqrt64Floor(uint64 value)
  26. {
  27.     //bshift = 8*8/2-1=32-1=31
  28.     uint64 bshift = ( CHAR_BIT * sizeof(bshift) >> 1 ) -1ULL;
  29.     uint64 tmp, ans = 0ULL;
  30.     size_t b=1U<<bshift;
  31.     do{
  32.         tmp = (((ans<<1)+b)<<bshift--);
  33.         if( value >= tmp) ans+=b, value-=tmp;
  34.     }while(b>>=1);
  35.     return ans;
  36. }
  37.  
  38. uint64 Sqrt64Ceil(uint64 value)
  39. {
  40.     //bshift = 8*8/2-1=32-1=31
  41.     uint64 bshift = ( CHAR_BIT * sizeof(bshift) >> 1 ) -1ULL;
  42.     uint64 tmp, ans = 0ULL, v=value;
  43.     size_t b=1U<<bshift;
  44.     do{
  45.         tmp = (((ans<<1)+b)<<bshift--);
  46.         if( v >= tmp) ans+=b, v-=tmp;
  47.     }while(b>>=1);
  48.     return ans += (ans * ans!=value);
  49. }
 
 
 
 
這篇就先到此，以後有 note 再於別篇補上。
arrow
arrow
    全站熱搜
    創作者介紹
    創作者 edisonx 的頭像
    edisonx

    Edison.X. Blog

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

    E-mail轉寄