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

使用bitset來取代bool陣列

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

const int N = 20000000; int sieve[(N>>5) + 1]; // 加一，以包含餘數。 int GET(int x) { return sieve[x>>5]>>(x&31) & 1; } void SET(int x) { return 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的倍數 int GET(int x) { return sieve[x>>5]>>(x&31) & 1; } void SET(int x) { return sieve[x>>5] |= 1<<(x&31); } int N2I(int i) { return (i-1)>>1; } int I2N(int i) { return (i<<1)+1; } #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); }

時間複雜度

考慮迴圈索引值i與j共有多少種：i一共有N種；j一共有(N/p1 - p1) + (N/p2 - p2) + ... + (N/sqrtN - sqrtN)種，其中p1 p2 ...是質數2 3 ...。時間複雜度就是兩者相加，為O(NloglogsqrtN) = O(NloglogN)。

1/1  + 1/2  + ... + 1/N     = O(logN)
1/p1 + 1/p2 + ... + 1/N     = O(loglogN)
1/p1 + 1/p2 + ... + 1/sqrtN = O(loglogsqrtN) = O(loglogN)

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

至於迴圈最內部，每個合數都被刪掉一次，總共O(N)次。

線性時間篩法

一邊製作質數表，一邊刪掉每個數的質數倍，時間複雜度就能達到O(N)。然而實際執行速度並沒有更快。

const int N = 20000000; bool sieve[N]; void eratosthenes() { 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; } } }

6n±1 Method

這是一個精簡版的篩法，原理是：只拿2和3這兩個質數先篩過一遍，剩下的數字則用試除法驗證是不是質數。

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+5也可以寫成6n-1，意義不變。）

6n-1到6n+1，再到下一個6n-1，再到6n+1，把這些要驗證的數字由小排到大，可以發現之間的差值會是2 4 2 4 2 4 ...不斷跳二跳四。實作程式碼時，就可以直接用加法加二加四， 而不必用乘法及加減法計算6n-1、6n+1，如此一來程式的執行效率會好一點。

驗證的順序是：數字2和3明顯的是質數，不必驗證；若是從數字5開始驗證，那麼下一個要驗證的數字就是5+2，再下一個就是5+2+4，再下一個就是5+2+4+2，如此不斷下去。

// 以試除法判斷質數 bool isprime(int n) { for (int i=5; i*i<=n; i+=2) if (n%i == 0) return false; return true; } // 只檢查6n±1。這些數字的間隔為2 4 2 4.... void make_prime() { cout << "找到質數2"; cout << "找到質數3"; for (int i=5, gap=2; i<1000000; i+=gap, gap=6-gap) if (isprime(i)) cout << "找到質數" << i; } vector<int> prime; // 以試除法判斷質數 bool isprime(int n) { for (int i=0; prime[i]*prime[i]<=n; i++) if (n % prime[i] == 0) return false; return true; } // 只檢查6n±1。這些數字的間隔為2 4 2 4.... void make_prime() { cout << "找到質數2"; prime.push_back(2); cout << "找到質數3"; prime.push_back(3); for (int i=5, gap=2; i<1000000; i+=gap, gap=6-gap) if (isprime(i)) {cout << "找到質數" << i; prime.push_back(i);} }

這個方法的時間複雜度是O(NsqrtN)，空間複雜度是O(1)。事實上6n±1法比篩法慢上許多。不過6n±1法的程式碼構造較為單純，當要枚舉的質數範圍不大時，有機會跑得比篩法還快。另外，6n±1法不需要開一條超大陣列來做計算，比起篩法它節省了很多空間，這也是它的優點。

Primality Test

質數測試，測試一個數字是否為質數。質數測試的方法相當多，其中有保證結果一定正確的方法，亦有結果不一定正確的方法。以下所介紹的，除了Divisibility Test的結果一定是正確的，其他方法的結果都不一定正確。

質數測試屬於P問題，不過以下介紹的皆非多項式時間的演算法，甚至是不保證結果正確的演算法。若對多項式時間、保證結果正確的演算法有興趣，可上網搜尋AKS Algorithm。

要進行質數測試，也可以直接用篩法把所有質數生出來，再來判斷質數。

Divisibility Primality Test

整除性測試法。依照質數定義，一個質數p不會被大於1且小於p的數字整除，只要把這些數字都拿來試除，就可以判定一個數字是不是質數。

撇開定義，這演算法其實就是窮舉所有可能的因數一一試除。

bool divisibiity_test(int n) { // 窮舉n所有可能的因數一一試除。 for (int d=2; d<n; ++d) if (n % d == 0) return false; // 不是質數 return true; // 是質數 } bool divisibiity_test(int n) { // 一個數字n不會有大於sqrt(n)的質因數（除了n本身以外） int sqrt_n = sqrt(n); for (int d=2; d<=sqrt_n; ++d) if (n % d == 0) return false; return true; }

這個演算法會進行sqrt(n)-1次除法，可推得時間複雜度為O(sqrt(n))，然而前提是：不管n多大，每次除法都是O(1)。

當要測試的數字很多時，可先建立質數表，進行質數測試時僅檢查質因數。

Fermat's Primality Test

費馬小定理：
若n是質數，則a^(n-1) ≡ 1 (mod n)。

費瑪質數測試法是運用費瑪小定理而想出的方法：

若n是質數，則費瑪小定理一定成立：a^(n-1) % n = 1一定成立。
若n是合數，則費瑪小定理不一定成立：a^(n-1) % n = 1有可能會成立。

當a^(n-1) % n = 1成立的時候就說n是質數。

這個演算法的結果不一定正確，有些合數也許會被判定成質數。但只要使用各式各樣的a來進行測試，那麼結果就會更加準確。

bool fermat_primality_test(int n, int a) { // 計算a^(n-1) % n的值，存於r中。 int r = 1; for (int i=1; i<n; ++i) r = r * a % n; return r == 1; }

Square Root Primality Test

若n是質數，則1~n之間，只有1與n-1，其平方後模n會餘1。
若ｎ是合數，則可能還會有其它的數字，其平方後模n會餘1。

（若以質數n為模，1的平方根只會等於±1，不會等於其他數。
　若以合數n為模，1的平方根只會等於±1，還可能會等於其他數。）

如果2~n-1的數字，其平方後模n都不餘1，就說n是質數。

這個演算法的結果不一定正確，有些合數也許會被判定成質數，例如22就會判定成質數，2^2、3^2、…、21^2模22後剛好都餘1。

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

Miller-Rabin Primality Test

這是結合Fermat's Primality Test與Square Root Primality Test的一個方法。這個演算法的結果不一定正確，有些合數也許會被判定成質數。

1. 選定一個底數a，用來進行費馬測試。
2. 把n-1拆解成m * 2^k的形式，m為奇數。
3. 觀察a^m模n。若為±1，
   則表示n最後可以通過費馬測試，接下來也不再出現平方根測試的反例。
   推定n為質數。
4. 依序觀察a^m、(a^m)^2 = a^(m * 2^1)、((a^m)^2)^2 = a^(m * 2^2)
   、…、((((a^m)^2)^2)…^2) = a^(m * 2^(k-1))這些數字：
　　甲、一旦發現有個數字平方後模n餘+1，
　　　　則表示n無法通過平方根測試。
　　　　推定n為合數。
　　乙、一旦發現有個數字平方後模n餘-1，
　　　　則表示n最後可以通過費馬測試，接下來也不再出現平方根測試的反例。
　　　　推定n為質數。
5. 最後無法判定結果。推定n為合數。

bool Miller_Rabin_primality_test(int n, int a) { // n-1 = m * 2^k int m = n – 1, k = 0; for (; !(m&1); m>>=1) k++; // t = a^m % n int t = pow(a, m) % n; if (t == +1 || t == -1) return true; for (int i=1; i<k; ++i) { t = t * t % n; if (t == +1) return false; if (t == -1) return true; } return false; }

步驟零。當t為±1：甲、t接下來的值會一直等於+1，於最後一步驟算式中可發現n通過費馬測試。乙、t接下來的值會一直等於1，我們不會檢查到違反平方根測試的情況。由甲乙兩點，推定n為質數。

步驟一到步驟k-1。一、當t為+1：n會通過費馬測試法，原理同前；但是n不會通過平方根測試法，因為上一個步驟的t不會是±1，而它的平方再模n（也就是此步驟的t）竟然等於1，也就是說我們發現了平方根測試的反例，故n必是合數。二、當t為-1：同步驟零的甲乙兩點，推定n是質數。三、若t為其他值：則無法判定n是質數或合數，繼續執行演算法。

步驟k。仍無法判定n是質數或是合數，故推定n為合數。

Factorization

因式分解。這裡談的是把一個正整數分解成質因數的連乘積：

n = 2^m1 * 3^m2 * 5^m3 * 7^m4 * 11^m5 * …

因數分解問題屬於NP問題，目前尚未有好解法。

Fundamental Theorem of Arithmetic（算數基本定理）

凡是正整數都可以藉由質因數分解成為一個獨一無二的式子，不同的n就會對應不同的(m1, m2, m3, …)，反方向亦同。

根據算數基本定理，凡是牽扯到乘法、因數、倍數的數學運算，都可以改變成比較簡單的運算。

分解前　　 | 分解後
-----------+--------------------------------------
乘除法     | 加減法
x * y      | (a1+b1, a2+b2, ...)
           |
整除       | 大於等於
x % y = 0  | (a1, a2, ...) ≥ (b1, b2, ...)
y | x      | a1≥b1 and a2≥b2 and ...
           |
最大公因數 | 最小值
gcd(x, y)  | (min(a1, b1), min(a2, b2), ...)
           |
最小公倍數 | 最大值
lcm(x, y)  | (max(a1, b1), max(a2, b2), ...)
           |
互質       | 對應項至少有一項為零
x⊥y       | min(a1, b1) = 0 and min(a2, b2) = 0 and ...
           | a1*b1 = 0 and a2*b2 = 0 and ...

算術基本定理闡述了另一種世界觀，把數字看作是質數的結合。質數的英文prime有著「原始就有」的意思，便是指質數是所有數字的根本。

Trial Division Factorization Method

把所有可能的因數拿來試除。用質因數會更好。

void trial_division(int n) { // 窮舉n所有可能的因數一一試除。 // 最簡單的方式是找由2到n當中的數字。 for (int d=2; d<=n; ++d) while (n % d == 0) { n / = d; cout << d; // 印出質因數 } } void trial_division(int n) { // 一個數字n不會有大於sqrt(n)的質因數（除了n本身以外） int sqrt_n = sqrt(n); for (int d=2; d<=sqrt_n; ++d) while (n % d == 0) { n / = d; cout << d; // 印出質因數 } if (n > 1) cout << n; // n是質數 }

當要因式分解的數字有很多筆時，可以先建好質數表，然後只拿質因數來試除。

UVa 516 583 10179 10290 10329 10392 10622 10780 10791

Fermat's Factorization Method

把一個數字分解成兩個數的乘積（這兩個數不一定會是質數）。一直遞迴分解下去，無法再分解的時候就找到質因數了。分解手法如下：

現在要找出 a 和 b 讓 n = a*b
使用公式 x^2 – y^2 = (x+y) * (x-y)
令 n = x^2 – y^2, a = x+y, b = x-y
窮舉整數 x，看看 sqrt(x^2-n) 是不是剛好就是一個整數 y，
如果是整數就找到一組 a 和 b 了。

void Fermat(int n) { for (int x = ceil(sqrt(n)); x<n; ++x) { int yy = x * x – n, y = sqrt(yy); if (y * y == yy) // 如果sqrt(x^2-n)是整數 { int a = x + y, b = x – y; cout << "n=" << a << "*" << b; // 遞迴分解a和b，就能完成因數分解。 Fermat(a); Fermat(b); } } }

Pollard's p-1 Factorization Method

這個方法可以找出n的其中一個質因數p。數學家Pollard證明當p-1的質因數都不大於b時，則p = gcd(2^b! - 1, n)。各位可以從wiki找到證明。:)

選定一個適當的b後，就可以進行演算法。有可能會找不到質因數。

void Pollard_p_minus_1(int n, int b) { int a = 2; // 計算2^b! - 1，存入a中。 for (int e=2; e <= b; ++e) e = expmod(a, e, n); // e = a^e % n p = gcd(a-1, n); // gcd(2^b! - 1, n) if (p == 1 || p == n) return -1; // 失敗 return p; }

Pollard's ρ Factorization Method

這個方法可以找出n的其中一個質因數p。有可能會失敗。

若有兩數x和y，p可以整除x-y，n不能整除x-y，則p = gcd(x-y, n)。

這個p可能是1、n、n的質因數，但是只要努力不懈地枚舉x和y，就有機會求得n的質因數。

以亂數產生器枚舉x和y

首先使用一個簡單的亂數產生器f(a_i+1) = a_i² + c，以及亂數種子a₀，其中c為整數常數。然後以(a₁ , a₀) (a₂ , a₁) (a₃ , a₂) … (a_n-1 , a_n-2)這些數對作為x和y。

int f(int a) { static int c = 1; // 不可以為0和-2 return a * a + c; } int PollardRho(int n) { static int a0 = 2; // 至少為+2 int x = a0, y = a0, p = 1; while (p == 1) { y = x; x = f(x); p = gcd(x-y, n); } if (p == n) return -1; // 失敗 return p; }

讀者可在此思考一下：

一、為何亂數產生器不選用f(a_i+1) = a_i + c這條更加簡單的式子？

二、為何a₀至少要是+2？

三、為何c不可以是0，也不可以是-2？讀者可將亂數產生器的式子，代入到x和y之中，計算一下x-y，並計算gcd(abs(x-y), n)。

四、如果我們將x和y對調，對結果有影響嗎？

五、為何p等於n的時候，就要結束迴圈，宣告失敗，不繼續找了？

小改進：x和y可以先模n，防止溢位

計算gcd(abs(x-y), n)的時候，其實就會拿abs(x-y)去模n。所以我們可以讓x和y每次都先模n，或者是讓亂數產生器每次都模n，這麼做不影響結果。不斷模n，可以防止x和y溢位，好處多多。

int PollardRho(int n) { static int a0 = 2; int x = a0, y = a0, p = 1; while (p == 1) { y = x; x = f(x) % n; // 快樂的模n，可防止溢位 p = gcd(abs(x-y), n); // 如果x比y小，就對調x和y } if (p == n) return -1; return p; }

小改進：偵測循環

方才的枚舉方式並不完美。我們知道，在模n的情況下，x和y最終必會重複出現，產生循環。我們得在循環開始時，就立刻結束演算法，否則接下來就會一直遇到重複的數字，進行重覆的運算。況且越早進入循環，就有越多數字白算。

Pollard採用了Floyd's Cycle Finding Algorithm的概念來偵測循環，以(a₁ , a₂) (a₂ , a₄) (a₃ , a₆) … (a_n , a_2n)這些數對作為x和y，一旦a_n等於a_2n即發現循環（當有模n的情況下），馬上結束演算法。

int f(int a) { static int c = 1; // c不可以為0或-2 return a * a + c; } int PollardRho(int n) { static int a0 = 2; // 至少為+2 int x = a0, y = a0, p = 1; while (p == 1) { x = f(x) % n; y = f(f(y) % n) % n; if (x == y) return -1; // 出現循環，馬上結束，宣告失敗 p = gcd(abs(x-y), n); } return p; }

演算法到此告一段落！各位在實作時，如果演算法失敗，可以更換a₀和c，或許就可以找出質因數了！

至於這個演算法的名稱由來，是因為不斷枚舉a₀ , a₁ , a₂ , ...，最終必會循環，繪圖時可以畫成一個ρ (rho)的形狀，因而得名。

非常可惜，沒有練習題！

Euler's Totient Function（Euler's φ Function）

這是一個公式。計算1到n的正整數當中，跟n互質（最大公因數是一）的數，總共有幾個。

首先要將n做質因數分解：

n = p1a1 × p2a2 × … × pkak   where p1 … pk are primes

以質因數計算Euler's Totient Function，φ唸做phi：

φ(n) = n × (1 - 1/p1) × (1 - 1/p2) × … × (1 - 1/pk)

可以採用這個順序計算，避免溢位：

n ÷ p1 × (p1-1) ÷ p2 × (p2-1) ÷ … ÷ pk × (pk-1)

如此就不用一個一個去計算最大公因數了，非常有效率！

當質因數分解採用試除法，計算Euler's Totient Function的時間複雜度為O(sqrtN)。預先建立質數表，得加速至O(π(sqrtN))，其中π(N)是1到N的質數個數。

UVa 10299 10179 11064

下面是Euler's Totient Function的一些性質，證明省略之：

φ(p) = p - 1                  iff p is prime.
φ(pa) = pa - pa-1              iff p is prime.
φ(n × m) = φ(n) × φ(m)        iff n and m are relatively prime.
φ(n) = φ(p1a1) × … × φ(pkak)   iff n = p1a1 × … × pkak
                              where p1 … pk are prime.

建立表格

未經改良的篩法，能求出每個數的質因數。運用篩法計算Euler's Totient Function，時間複雜度為O(NloglogN)。

int phi[10000]; void phi_table() { phi[1] = 1; for (int i=2; i<10000; ++i) if (!phi[i]) for (int j=i; j<10000; j+=i) { if (!phi[j]) phi[j] = j; phi[j] = phi[j] / i * (i-1); } }

或者，首先運用篩法，為每個數求得一個質因數；然後運用Euler's Totient Function的性質，實施Dynamic Programming。此時得以套用線性時間篩法，時間複雜度為O(N)。

int phi[10000]; void phi_table() { // 篩法。為每個數求得一個質因數。O(NloglogN)。 // 亦可運用線性時間篩法。O(N)。 for (int i=2; i<100; ++i) if (!phi[i]) for (int k=(10000-1)/i, j=i*k; k>=i; k--, j-=i) if (!phi[k]) phi[j] = i; // i是j的其中一個質因數 // 計算Euler's Totient Function。O(N)。 phi[1] = 1; for (int i=2; i<10000; ++i) if (!phi[i]) // i是質數 phi[i] = i-1; else // i不是質數 { int j = i / phi[i]; // i分解成j*phi[i] if (j % phi[i] == 0)// j與phi[i]不互質 phi[i] = phi[j] * phi[i]; else phi[i] = phi[j] * (phi[i] - 1); } }

UVa 10820 10990 11327 11424 11426

延伸閱讀：最大公因數的傅立葉轉換就是Euler's Totient Function

http://en.wikipedia.org/wiki/Euler%27s_totient_function#Fourier_transform

       N-1
φ(n) =  Σ { gcd(n, k) ÷ ei*2π*(n/N)*k }
       k=0

Möbius Function

用排容原理求一個數字的所有因數總和。

ICPC 2116