使用加法，徹底分解數字

例如5 = 1 + 1 + 1 + 1 + 1。

不可分解的單元是1。不稀奇。

使用乘法，徹底分解數字

加法換成乘法之後，事情變化多端，例如5 = 5，例如6 = 2 × 3。可以發現，2、3、5、7、11等，是不可分解的單元。至於4、6、8、9、10等，可以再分解。

不可分解的單元叫做「質數」！非常稀奇！

接下來要介紹的演算法有：從小到大列出質數（建立質數表）、判斷一個數是不是質數（質數測試）、使用乘法湊得給定的數（質因數分解）。

prime generation

製作質數表。

  2   3   5   7  11  13  17  19  23  29
 31  37  41  43  47  53  59  61  67  71
 73  79  83  89  97 101 103 107 109 113
127 131 137 139 149 151 157 163 167 173
179 181 191 193 197 199 211 223 227 229
233 239 241 251 257 263 269 271 277 281
......

sieve of Eratosthenes

這是一個製作質數表的演算法。簡稱「篩法」。

列出所有正整數。從2開始，刪掉2的倍數。找下一個未被刪掉的數字，找到3，刪掉3的倍數。找下一個未被刪掉的數字，找到5，刪掉5的倍數。如此不斷下去，就能刪掉所有合數，找到所有質數。

質數有無限多個，證明省略。我們無法找到所有質數，通常是預先訂立一個範圍，只找到範圍內的所有質數。

bool prime[20000000]; void eratosthenes() { for (int i=0; i<20000000; i++) // 初始化 prime[i] = true; prime[0] = false; // 0 和 1 不是質數 prime[1] = false; // 找下一個未被刪掉的數字 for (int i=2; i<20000000; i++) if (prime[i]) // 刪掉質數i的倍數，從兩倍開始。保留原本質數。 for (int j=i+i; j<20000000; j+=i) prime[j] = false; }

欲刪掉質數i的倍數之時，早已刪掉1倍到i-1倍了，直接從i倍開始。

bool prime[20000000]; void eratosthenes() { for (int i=0; i<20000000; i++) prime[i] = true; prime[0] = false; prime[1] = false; for (int i=2; i<20000000; i++) if (prime[i]) // 刪掉質數i的倍數，從i倍開始。小心溢位。 for (int j=i*i; j<20000000; j+=i) prime[j] = false; }

欲刪掉質數i的倍數之時，早已刪掉「小於i的質數、其倍數」倍了，直接刪掉「大於等於i的質數、其倍數」倍。

乍看下程式碼增多而變慢，實際上cache miss減少而變快。

bool prime[20000000]; void eratosthenes() { for (int i=0; i<20000000; i++) prime[i] = true; prime[0] = false; prime[1] = false; for (int i=2; i<20000000; i++) if (prime[i]) // k是倍率，j是質數i的倍數。 // 由大到小判斷，當prime[k]成立時， // k才能涵蓋所有「大於等於i的質數、其倍數」倍。 for (int k=(20000000-1)/i, j=i*k; k>=i; k--, j-=i) if (prime[k]) prime[j] = false; } // 運用小範圍的prime[k]，避免存取大範圍的prime[j]， // 大幅減少cache miss。

一個合數x，必定有一個小於等於sqrt(x)的質因數。所有小於等於sqrt(x)的質數，刪掉這些質數的倍數，就能刪掉所有合數了。

bool prime[20000000]; void eratosthenes() { for (int i=0; i<20000000; i++) prime[i] = true; prime[0] = false; prime[1] = false; // 只需要刪掉sqrt(20000000)以下的質數的倍數。 int sqrt_20000000 = sqrt(20000000); for (int i=2; i<=sqrt_20000000; i++) if (prime[i]) for (int k=(20000000-1)/i, j=i*k; k>=i; k--, j-=i) if (prime[k]) prime[j] = false; }

顛倒true和false，節省初始化時間。

// 全域變數的陣列，每一格將自動初始化為false。 bool sieve[20000000]; void eratosthenes() { int sqrt_20000000 = sqrt(20000000); sieve[0] = sieve[1] = true; for (int i=2; i<=sqrt_20000000; i++) if (!sieve[i]) for (int k=(20000000-1)/i, j=i*k; k>=i; k--, j-=i) if (!sieve[k]) sieve[j] = true; }

製作質數表。篩法結束之後，掃描一次陣列即可。

bool sieve[20000000]; void eratosthenes() { ......; // 製作質數表 vector<int> prime; for (int i=2; i<20000000; i++) if (!sieve[i]) prime.push_back(i); }

UVa 406 516 524 543 10140 10311 10394 11408

使用bitset來取代bool陣列

一個int有32個位元，可以當作32個欄位來使用，節省記憶體空間，減少cache miss。

const int N = 20000000; int sieve[(N>>5) + 1]; // 加一，以包含餘數。 // inline關鍵字讓編譯器用「函式內容」覆蓋掉「函式呼叫」。 inline int GET(int x) { return sieve[x>>5]>>(x&31) & 1; } inline void SET(int x) { sieve[x>>5] |= 1<<(x&31); } void eratosthenes() { // register關鍵字讓變數存放於CPU register，加快存取速度。 register int i, j, k; int sqrt_N = sqrt(N); SET(0); SET(1); for (i=2; i<=sqrt_N; i++) if (!GET(i)) // 因為陣列已經加一，N-1改成N也無妨。 for (k=N/i, j=i*k; k>=i; k--, j-=i) if (!GET(k)) SET(j); }

不處理2的倍數

不處理2的倍數，節省一半記憶體，增進一點速度。

令陣列第0格代表數字1、陣列第1格代表數字3、陣列第2格代表數字5、……，以此類推。

const int N = 1<<25; int sieve[(N>>6) + 1]; // 不處理2的倍數 // C++風格 inline int GET(int x) { return sieve[x>>5]>>(x&31) & 1; } inline void SET(int x) { sieve[x>>5] |= 1<<(x&31); } inline int N2I(int i) { return (i-1)>>1; } inline int I2N(int i) { return (i<<1)+1; } // C風格 #define GET(x) (sieve[(x)>>5]>>((x)&31) & 1) #define SET(x) (sieve[(x)>>5] |= 1<<((x)&31)) #define N2I(i) (((i)-1)>>1) // 索引值 = (數字 - 1) / 2 #define I2N(i) (((i)<<1)+1) // 數字 = 索引值 * 2 + 1 // 驗證一個數字是不是質數 bool isprime(int x) { return x==2 || x>2 && (x&1) && !GET(N2I(x)); } void eratosthenes() { int half_sqrt_N = sqrt(N) / 2; int half_N = N >> 1; register int i, j, k; // 從3開始。N2I(3) = 1。 for (i=1; i<=half_sqrt_N; i++) if (!GET(i)) for (j=2*i*(i+1), k=2*i+1; j<half_N; j+=k) SET(j); } void eratosthenes() { int half_sqrt_N = sqrt(N) / 2; register int i, j, k, ii; for (i=1; i<=half_sqrt_N; i++) if (!GET(i)) for (ii=I2N(i), j=N2I(N/ii), k=ii*j+i; j>=i; j--, k-=ii) if (!GET(j)) SET(k); }

時間複雜度

考慮內層迴圈索引值j一共有多少種：(N/2 - 2) + (N/3 - 3) + (N/5 - 5) + ... + (N/sqrtN - sqrtN) = O(NloglogN)。

1/1 + 1/2 + 1/3 + ... + 1/N     = O(logN)
1/2 + 1/3 + 1/5 + ... + 1/N     = O(loglogN)
1/2 + 1/3 + 1/5 + ... + 1/sqrtN = O(loglogsqrtN) = O(loglogN)

https://en.wikipedia.org/wiki/Divergence_of_the_sum_of_the_reciprocals_of_the_primes
1/1 + 1/2 + 1/3 + ... + 1/n - ln(n)     趨近於Euler–Mascheroni constant
1/2 + 1/3 + 1/5 + ... + 1/n - ln(ln(n)) 趨近於Meissel–Mertens constant

segmented sieve of Eratosthenes

篩法當中，以不大於sqrtN的質數，篩選大於sqrtN的數。

篩法需要大量記憶體空間。為了節省記憶體空間，於是有了分段處理的手法：一、首先求得不大於sqrtN的質數。二、將記憶體切成小段，每段長度sqrtN。分別篩選每一段。

空間複雜度從O(N)變成O(sqrtN)，大幅減少cache miss。

http://www.geeksforgeeks.org/segmented-sieve/

也許可以作為平行演算法的經典範例。

一次方時間的篩法

一邊製作質數表，一邊刪掉每個數的質數倍。每個合數只遇到最小質因數（次小因數），每個合數只讀取一次。時間複雜度O(N)。

const int N = 20000000; bool sieve[N]; void linear_sieve() { vector<int> prime; for (int i=2; i<N; i++) { if (!sieve[i]) prime.push_back(i); for (int j=0; i*prime[j]<N; j++) { sieve[i*prime[j]] = true; if (i % prime[j] == 0) break; } } } const int N = 20000000; bool sieve[N]; void linear_sieve() { vector<int> prime; for (int i=2; i<N; i++) { if (!sieve[i]) prime.push_back(i); for (int p: prime) { if (i * p >= N) break; sieve[i * p] = true; if (i % p == 0) break; } } }

primality test

「質數測試」。測試一個數字是否為質數。

質數測試屬於P問題。然而以下介紹的都是指數時間演算法，甚至是不保證結果正確的演算法。

多項式時間、保證結果正確的演算法，請逕行上網搜尋「AKS primality test」；但由於方法繁雜，實務上沒人使用。

質數測試，可以直接以篩法建立質數表，再來判斷質數。然而建立質數表需要大量記憶體，因此又發明了其他方法。

divisibility primality test

整除性測試法即是試除法。

質數無法以乘法分解，質數沒有任何因數（除了1和自身）。窮舉所有可能的因數，一一試除。時間複雜度O(n)。

bool isprime(int n) { if (n <= 1) return false; // 窮舉所有可能的因數，一一試除（除了1和自身） for (int d=2; d<n; ++d) if (n % d == 0) return false; return true; }

因數成雙成對（平方根跟自己一對）。檢查了小於等於平方根的因數，就不必再檢查大於平方根的因數。時間複雜度O(sqrt(n))。

bool isprime(int n) { if (n <= 1) return false; // 窮舉小於等於sqrt(n)的因數，一一試除（除了1） for (int d=2; d*d<=n; ++d) if (n % d == 0) return false; return true; }

2n+1 method

因數可以分解成質數。取質數試除，效果等同於取因數試除。反過來說，大可不必取合數試除。

2的更大倍數是合數，不必試除，節省一半時間。

bool isprime(int n) { if (n <= 1) return false; if (n == 2) return true; // 首先以2試除 if (n % 2 == 0) return false; // 2的更大倍數是合數，不必試除。 for (int d=3; d*d<=n; d+=2) // 跳二 while (n % d == 0) return false; return true; }

6n±1 method

2和3的最小公倍數是6。所有數字分為6n、6n+1、6n+2、6n+2、6n+3、6n+4、6n+5六種（n是倍率）。除以六餘零的數字為6n，除以六餘一的數字為6n+1，以此類推。

6n、6n+2、6n+3、6n+4是2或3的倍數，不是質數。只需要拿6n+1和6n+5試除。

從6n+1，跳到6n+5，跳到下一個6n+1，再跳到6n+5，差值是4 2 4 2 4 2...，不斷跳四跳二。

6n+5改寫成6n-1，意義不變。6n+1與6n-1簡寫成6n±1。

2和3預先處理。從5開始試除，下一個是5+2，再下一個是5+2+4，再下一個是5+2+4+2，以此類推。

實作程式碼，用加法加二加四，不用乘法計算6n-1、6n+1，效率較佳。

bool isprime(int n) { if (n <= 1) return false; if (n == 2) return true; if (n == 3) return true; // 首先以2和3試除 if (n % 2 == 0) return false; if (n % 3 == 0) return false; // 然後只檢查6n±1。這些數字的間隔為2 4 2 4... for (int d=5, gap=2; d*d<=n; d+=gap, gap=6-gap) if (n % d == 0) return false; return true; }

wheel factorization

6n±1法進行推廣，使用頭幾個質數。

取2 3 5的最小公倍數30，建立30n+k法。取2 3 5 7的最小公倍數210，建立210n+k法。取所有質數，等同篩法。

取越多質數，即可剔除越多合數，事前節省時間。

然而，取越多質數，也會建立越多規則，事中浪費時間。

這是取捨兩難的情況，兩者之間必須取得平衡。一般來說，大家只取2 3或者2 3 5。

更有甚者，當你還在煩惱要取多少質數，普通的試除法早就已經執行完畢了。更有甚者，當你還在煩惱要如何實作試除法，網路上已經有現成的質數表了。所謂的取得平衡，端看你的視野高度。

square root primality test

餘數系統的平方根有個特性：

以質數n為模數，0到n-1之間，只有1與n-1，平方之後等於1。
以質數n為模數，1的平方根一定是1與n-1（即±1）。

平方根測試法，運用了此特性：

以n為模數，當2到n-2的每一個數字，平方之後皆不等於1，就推定n是質數。

此演算法的結果不一定正確。通過測試的數字，可能是合數或質數；無法通過測試的數字，一定是合數。

以合數n為模數，0到n-1之間，除了1與n-1，還可能有其它數字，平方之後等於1。
以合數n為模數，1的平方根可能不只是1與n-1（不只±1）。

有些合數會被誤判成質數，例如22會被誤判成質數，2²、…、20²模22之後皆不等於1。

bool square_root_primality_test(int n) { for (int i=2; i<=n-2; ++i) if (i * i % n == 1) return false; return true; }

Fermat primality test

費瑪小定理：

若n是質數，a小於n，則aⁿ ≡ a (mod n)。
若n是質數，a小於n，a不是零，則aⁿ⁻¹ ≡ 1 (mod n)。

費瑪質數測試法，運用了費瑪小定理：

n是質數，費瑪小定理一定成立：aⁿ⁻¹ % n = 1一定成立。
n是合數，費瑪小定理可能成立：aⁿ⁻¹ % n = 1可能成立。

當aⁿ⁻¹ % n = 1成立，就推定n是質數。

此演算法的結果不一定正確。通過測試的數字，可能是合數或質數；無法通過測試的數字，一定是合數。

使用各式各樣的a來實施測試，那麼判定結果就更準確。

bool fermat_primality_test(int n) { // 隨便取得一個數值當作a， // 但是必須大於1、小於n，才有意義。 // srand(time(0)); int a = rand() % (n-2) + 2; // 計算aⁿ⁻¹ % n的值，存於x中。 int x = 1; for (int i=1; i<n; ++i) x = x * a % n; return x == 1; // return pow(a, n-1, n) == 1; }

Miller–Rabin primality test

結合square root primality test與Fermat primality test。

此演算法的結果不一定正確。通過測試的數字，可能是合數或質數；無法通過測試的數字，一定是合數。

一、選定一個底數a，大於1、小於n，用來進行費瑪測試。
　　待測數字n，會是費瑪測試的模數。
二、令 n-1 = u × 2ᵗ，其中t盡量大（u為奇數）。
三、觀察a^u這個數字：
　　若等於±1，
　　則表示步驟四，所有數字都是1，永不出現平方根測試的反例。
　　也表示步驟五，最終將通過費瑪測試。
　　推定n是質數。
四、依序觀察a^(u × 2¹)、a^(u × 2²)、…、a^(u × 2ᵗ⁻¹)這些數字，
　　每個數字正好是前一個數字的平方：
　甲、一旦發現有個數字的平方等於+1，
　　　則表示無法通過平方根測試。
　　　（但是步驟五，最終將通過費瑪測試。）
　　　確定n是合數。
　乙、一旦發現有個數字的平方等於-1，
　　　則表示接下來的數字都是1，永不出現平方根測試的反例。
　　　也表示步驟五，最終將通過費瑪測試。
　　　推定n是質數。
五、觀察a^(u × 2ᵗ)這個數字：
　甲、若等於+1，表示通過費瑪測試，推定n是質數。
　乙、若不等於+1，表示無法通過費瑪測試，確定n是合數。
　回、也就是說，a^(u × 2ᵗ⁻¹)必須等於±1，平方之後才會等於+1。
　　　步驟五才能通過費瑪測試。
　回、由於步驟四已經檢查過a^(u × 2ᵗ⁻¹)是否為±1，
　　　所以步驟五可以省略，直接確定n是合數。

bool miller_rabin(int n) { if (n < 2) return false; // srand(time(0)); int a = rand() % (n-2) + 2; int u = n-1, t = 0; while (u % 2 == 0) u >>= 1, t++; int x = pow(a, u, n); // x = a ^ u % n; if (x == 1 || x == n-1) return true; for (int i=0; i<t-1; i++) { x = mul(x, x, n); // x = x * x % n; if (x == 1) return false; if (x == n-1) return true; } return false; }

步驟二：不斷分解因數2，最多log(n)步，需時O(log(n))。

步驟三：餘數次方，divide-and-conquer method，需時O(log(n))。

步驟四：根據步驟二，最多log(n)-1個數字，需時O(log(n))。

總時間複雜度O(log(n))。

strong probable-prime base（SPRP）

以過濾合數的觀點來看，多取幾個相異的底數a實施測試，判定結果就更準確。

事實上已經有熱心人士，找出特定的底數組合，仔細檢查了各種數字的判定結果，保證判定結果正確。例如底數組合{2, 7, 61}就能正確判斷2³²以下的數字是不是質數。

bool isprime(int n) { // 預先判斷偶數與1，節省一點時間。 if (n == 2) return true; if (n < 2 || n % 2 == 0) return false; int u = n - 1, t = 0; while (u % 2 == 0) {u >>= 1; t++;} // 推定是質數，就實施下一次測試； // 確定是合數，就馬上結束。 int sprp[3] = {2, 7, 61}; // 足以涵蓋int範圍 for (int k=0; k<3; ++k) { // 預先判斷底數是否為簡易數字，節省一點時間。 int a = sprp[k] % n; if (a == 0 || a == 1 || a == n-1) continue; long long x = pow(a, u, n); if (x == 1 || x == n-1) continue; for (int i = 0; i < t-1; i++) { x = mul(x, x, n); if (x == 1) return false; if (x == n-1) break; } if (x == n-1) continue; return false; } // 通過全部測試，確定是質數。 return true; }

UVa 10956

Lucas primality test

Baillie–PSW primality test

目前仍然可以正確判斷質數，尚未找到反例。