--- tags: 成大高階競技程式設計 2021 --- # 2021 Week14 - Number Theory ## 最大公因數(Greatest Common Divisor) 每個整數都會有一個以上的因數,而兩個整數可能會有一些共通的因數, 其中最大的因數就被稱為最大公因數。 雖然 STL 中有[現成的函式](https://hackmd.io/@XDEv11/NCKU-AdvCP-2021-CXX#stdgcd--stdlcm-C17)可以直接用,不過實作的演算法也很值得學習。 ### 基本性質 1. $\gcd(a, b) = \gcd(b, a)$ 2. $\gcd(a, b) = \gcd(a \ \% \ b, b)$ 3. $\gcd(a, 0) = | a |$ > $\%$ 為 C/C++ 中的取餘數運算 這邊只解釋第 $2$ 點,令 $g = \gcd(a, b)$,則 $a$ 和 $b$ 分別可以表示為 $a = a'g$ 以及 $b = b'g$ 模運算的算式就可以被改寫為 $a \ \% \ b = a - \lfloor \frac{a}{b} \rfloor b = (a' - \lfloor \frac{a}{b} \rfloor b')g$,可以發現 $a \ \% \ b$ 也是 $g$ 的倍數 並且因為 $a'$ 與 $b'$ 互質,所以 $a' - \lfloor \frac{a}{b} \rfloor b'$ 也會與 $b'$ 互質 從上述可得知 $\gcd(a, b) = \gcd(a \ \% \ b, b)$ ### 輾轉相除法(Euclidean Algorithm) 整合前面[基本性質](#基本性質)有提到的 $3$ 種性質而成的演算法 ```cpp! int gcd(int a, int b) { while (b) { a %= b; // 性質 2 swap(a, b); // 性質 1 } return abs(a); // 性質 3 } ``` 注意到因為每次 $b$ 都會比前一次的 $b$ 還要接近 $0$,所以不會形成無窮迴圈 > 這邊寫接近 $0$ 而不是小於前一次的 $b$ 是因為在 C/C++ 中如果被除數是負數的話 > 模運算得到的結果也會是負數 ### 貝祖定理(Bézout's Lemma) $ax + by = m$ 在上列等式中,如果 $m$ 是 $\gcd(a, b)$ 的倍數,則 $x, y$ 有無限多組**整數**解,否則為無整數解。 :::spoiler Proof 證明過程與[貝祖定理](https://zh.wikipedia.org/wiki/%E8%B2%9D%E7%A5%96%E7%AD%89%E5%BC%8F#%E6%95%B4%E6%95%B0%E4%B8%AD%E7%9A%84%E8%A3%B4%E8%9C%80%E5%AE%9A%E7%90%86)中寫到的過程一樣,這邊只是稍微講得比較詳細一點而已, * 如果 $a$ 為 $0$,$\gcd(a, b) = b$ 等式變成 $by = m$ 如果 $m$ 是 $b$ 的倍數的話 $x, y$ 就有無限多組整數解,因為 $x$ 可以是任意數, 否則為無整數解 > $b = 0$ 的情況同理 * 若 $a, b$ 都不為 $0$,令 $\gcd(a, b) = g$ 定義 $A = \left\{ ax + by \ | \ (x, y) \in \mathbb{Z}^2 \right\}$ 首先先來證明 $A$ 中最小的正整數 $d_0$ 就是 $g$ * $d_0 \leq g$ 因為 $d_0$ 存在於 $A$ 中,因此一定存在 $x_0, y_0$ 使得 $ax_0 + by_0 = d_0$ 成立 令 $P$ 為 $A$ 中任意的**正整數**,$P = ax + by$ $P$ 對 $d_0$ 做除法後可以得到一個商 $q$ 與一個餘數 $r$ 則 $P$ 就可以表示成 $P = qd_0 + r$ 將等式代換一下就會變成 $ax + by = q(ax_0 + by_0) + r$ 移項就可以得到 $a(x - qx_0) + b(y - qy_0) = r$ 這邊可以發現 $r$ 也是屬於 $A$ 的元素,但是要注意到 $r$ 是除法得到的餘數 因此必須要符合 $0 \leq r < d_0$,而 $d_0$ 又是 $A$ 中**最小**的正整數, 所以 $r$ 必為 $0$ $r = 0$ 就代表 $P$ 是 $d_0$ 的倍數,因此就可以得知 $A$ 中任意的正整數都是 $d_0$ 的倍數 $A$ 中所有的負整數也都是 $d_0$ 的倍數,因為這些負整數乘以 $-1$ 後就是 $A$ 中的某個正整數 而 $a$ 和 $b$ 也都是 $A$ 中的某個整數($(x, y)$ 分別代入 $(1, 0)$ 以及 $(0, 1)$) 因此 $d_0$ 是 $a$ 與 $b$ 的公因數,也可以得知 $d_0 \leq g$ * $d_0 \geq g$ 因為 $g$ 是 $a, b$ 的最大公因數,所以我們可以寫成 $a = a'g$ 以及 $b = b'g$ 而 $d_0$ 的定義就可以寫成 $d_0 = ax_0 + by_0 = (a'x_0 + b'y_0)g$ 也就是 $d_0$ 是 $g$ 的倍數,也代表 $d_0 \geq g$ 因為同時符合 $d_0 \leq g$ 以及 $d_0 \geq g$,因此可以得證 $d_0 = g$ 從上面證明過程中可以知道如果 $m$ 不為 $g$ 的倍數,代表 $m$ 不在 $A$ 之中,也就代表 $x, y$ 無整數解 反之如果 $m$ 是 $g$ 的倍數,代表 $m$ 在 $A$ 中,並且因為是二元一次多項式,所以有無限多組整數解 ::: ### 擴展歐幾里德算法(Extended Euclidean Algorithm) 解方程式 $ax + by = \gcd(a, b)$,除了算出 $\gcd(a, b)$ 之外 還找出一組 $(x, y)$ 使得上面等式成立 做法就是在[輾轉相除法](#輾轉相除法(Euclidean-Algorithm))的過程中,順便計算 $(x, y)$ ```cpp! int gcd(int a, int b, int &x, int &y) { if (b == 0) { x = a < 0 ? -1 : 1; y = 0; // 可以設為任意數 return x * a; } int g = gcd(b, a % b, y, x); y -= a / b * x; return g; } ``` $$ r = a \ \% \ b \implies r = a - \lfloor \frac{a}{b} \rfloor b $$ $$ \begin{array}{rl} by + rx &= by + (a - \lfloor \frac{a}{b} \rfloor b)x \\&= ax + by - \lfloor \frac{a}{b} \rfloor bx \\&= ax + b(y - \lfloor \frac{a}{b} \rfloor x) \end{array} $$ ## 同餘運算(Modular Arithmetic) 在競賽中若計算結果超出了數值範圍外,則通常會要求結果對某個數字**取餘數** 所以得熟悉取完餘數後的數字們之間的**關係**及**運算**。 若**數字們**都對 $n$ **取餘數**後,它們就處於"**模 $n$**"的狀態下了 將 $a$ 除以 $n$ 的餘數記做 $a \space (\text{mod } n)$ > 記得檢查若 a 為**負數**的情況 若 $a \space (\text{mod } n)$ 與 $b \space (\text{mod } n)$ 相等,則記 $a \equiv b \space (\text{mod } n)$,稱 $a, b$ 有**同餘關係** 也就是說 $\forall k \in \mathbb{Z}\centerdot a \equiv a + k\cdot n \space (\text{mod } n)$ > 可將括號內看作數字們處於一個**特定狀態**下,而模 $n$ 通常會省略括號 並且若 $$ a \equiv b \mod n \\c \equiv d \mod n $$ 則 $$ a+c \equiv b+d \mod n\\a\times c \equiv b\times d \mod n $$ ### 歐拉定理 如果 $a, n$ **互質**[^coprime],那麼 $$ a^{\phi(n)} \equiv 1 \mod n $$ 這裡 $\phi(n)$ 是與 $n$ 互質且小於 $n$ 的**正整數**的個數 > 例如 $\phi(6) = 2$ 因為 $1, 5$ 和 $6$ 互質 ### 費馬小定理 如果 $a, p$ 互質,且 $p$ 為**質數**,那麼 $$ a^{p-1} \equiv 1 \mod p $$ > 是歐拉定理的特例 ### 模反元素 數字們模 $n$ 後,對於加法、減法、乘法,其性質都與以往差不多,但對數字做**除法**卻不盡人意 > 在同餘運算中只會有**整數**,不會有浮點數 **反元素**指的是元素與其**運算**後為**單位元素**的元素 例如: - 整數 $a$ 的**加**法反元素為 $-a$ - 整數 $a$ 的**乘**法反元素為 $1\over a$ 而元素 $a$ 的**模 $n$ 反元素**以 $a^{-1}$ 表示,滿足 $$ a \cdot a^{-1} \equiv 1 \mod n $$ 根據 [Bézout’s Lemma](#貝祖定理(Bézout’s-Lemma)) - $ax + ny = \gcd(a, n) = 1$ 在模 $n$ 有 $ax \equiv 1 \mod n$ - $ax + ny = \gcd(a, n) \not= 1$ 在模 $n$ 有 $ax \not\equiv 1 \mod n$ 也就是說 $\gcd(a, n) = 1$ (互質),$a$ 才擁有模 $n$ 反元素 #### 求出模 $n$ 反元素 - 當 $n$ 為**質數**時: 根據[費馬小定理](#費馬小定理)有 $$ a \cdot a^{n-2} \equiv 1 \mod n $$ 所以 $a^{n-2}$ 是 $a$ 的模反元素 > 下面教的[**快速冪**](#快速冪(Exponentiating-by-Squaring))能在 $O(\log_2 n)$ 得到 $a^{n-2}$ - 當 $n$ **非質數**時: 根據 [Bézout’s Lemma](#貝祖定理(Bézout’s-Lemma)) 有 $$ \begin{split} &a\cdot x + n\cdot y \space &= 1 \\ \Rightarrow \space &a \cdot x &\equiv 1 \mod n \end{split} $$ 可用[擴展歐幾里得演算法](#擴展歐幾里德算法(Extended-Euclidean-Algorithm))找到 $x = a^{-1}$ ## 質數(Prime Numbers) 若數 $n > 1$ 能被 $1$ 與 $n$ 以外的數整除,則 $n$ 為合數 ### $[2, \sqrt{n}]$ 判斷法 假設 $n$ 是**合數**,則 $n = x\cdot y$ 其中 $1 < x \le y$,可證明 $x$ 的大小不超過 $\sqrt{n}$; 所以若有 $x \in [2, \sqrt{n}]$ 滿足 $x\mid n$,則 $n$ 不是質數。 > $x\mid n$ 表示 $x$ 整除 $n$ ```cpp for (int i = 2; i <= sqrt(n); i++) if (n%i == 0) return false; return true; ``` ### Fermat's Primality Test 根據[費馬小定理](#費馬小定理), $$ n \text{ is prime } \Rightarrow a^{n-1} \equiv 1 \mod n \text{, where }n \not\mid a $$ 雖然根據**邏輯規則**以下**逆命題不見得成立** $$ n \text{ is prime } \Leftarrow a^{n-1} \equiv 1 \mod n \text{, where }n \not\mid a $$ 但卻有很**高的機率**[^probability_a]是成立的! ```cpp int a = max(rand()%n, 2); // a in [2, n-1] int x = power(a, n-1, n); // 求 $a^{n-1}$ 除以 n 的餘數 return x == 1; ``` > 如果 `power` 函數是[快速冪演算法](#快速冪(Exponentiating-by-Squaring)),則複雜度 $O(\log n)$ 根據**邏輯規則**,如果回傳 `false`,那麼保證 $n$ 不是質數 > $p \Rightarrow q$ 與 $\neg p \Leftarrow \neg q$ 是等價的 注意這是**測試法**,有時會遇到難**判斷**的數,例如 $n = 561$ 對大部分的 $a$ 回傳 `true` ### Miller Rabin Primality Test 繼續從 [Fermat's Primality Test](#Fermat’s-Primality-Test) 發展: 除 $2$ 以外的質數都為奇數,故**只關注奇數** $n$ 就好 則 $n-1$ 是個**偶數**,可將其寫成 $n-1 = 2^t \cdot d$ 再注意到,若 $i < t$ 有 $a^{2^i \cdot d} \equiv \pm 1 \mod n$ 則 $a^{2^t \cdot d} \equiv 1 \mod n$ > 因為 $({\pm 1})^2$ 等於 $1$ 嘛 並且若 $n$ 是個質數,則 $$ \begin{split} x^2 \equiv 1 \mod n &\Rightarrow n \mid x^2 - 1 = (x+1)\cdot (x-1) \\ &\Rightarrow n \mid x+1 \lor n \mid x-1 \\ &\Rightarrow x \equiv -1 \mod n \lor x \equiv 1 \mod n \end{split} $$ > 如果 $n$ 是合數,第二列 $"\Rightarrow"$ **不一定**會成立 有了上述性質,又能使質數測試成功機率大大提升 ```cpp if (n%2 == 0) return n == 2; // 2 為質數,反之偶數非質數 a %= n; if (a == 0) return true; int d = n-1, t = 0; while (d%2 == 0) t++, d /= 2; int x = power(a, d, n); while (t--) { int nx = power(x, 2, n); if (nx == 1 && x != 1 && x != n-1) return false; x = nx; } return x == 1; ``` > 根據**邏輯規則**,如果回傳 `false`,那麼保證 $n$ 不是質數 經測試若數 $n$ 在下列範圍,皆能**正確判斷**質數: - $n < 2^{32}$ 代入所有 $a \in \{2, 7, 61 \}$ - $n < 2^{64}$ 代入所有 $a \in \{2, 325, 9375, 28178, 450775, 9780504, 1795265022 \}$ 這樣就能在 $O(\log n)$ 完美判斷常用質數了! ## 質因數分解 根據唯一分解定理,任何合數都能被分解成一些質數的積 > 算術基本定理 ### $[2, \sqrt{n}]$ 試除法 > 與[$[2, \sqrt{n}]$ 判斷法](#2-sqrtn-判斷法)原理相同 合數 $n = p_1^{t_1} \cdot p_2^{t_2} \cdot \cdots \cdot p_k^{t_k}$ 為多個質數的乘積 不失一般性的有 $p_1 < p_2 < \cdots < p_{k-1} < \sqrt{n}$,所以在範圍 $[2, \sqrt{n}]$ **從小至大**找數試除即可。 ```cpp for (int p = 2; p <= sqrt(n); p++) { int t = 0; while (n%p == 0) n /= p, t++; if (t) factors.push_back({p, t}); } if (n != 1) factors.push_back({n, 1}); ``` ## 質數篩檢 顧名思義,通常會想知道在一定的範圍內**哪些數是質數** ### Sieve of Eratosthenes 其精神是將質數的**倍數**都設為**非質數**(合數) ![](https://i.imgur.com/aws94bk.gif) ```cpp vector<bool> is_p(maxn, true); is_p[1] = false; for (int n = 2; n < sqrt(maxn); n++) { if (!is_p[n]) continue; for (int m = n*n; m < maxn; m+=n) is_p[m] = false; } ``` `n*n` 為對於所有 $x < n$,在數到 $n$ 之前 $x\cdot n$ 已經被篩過了。 >記得檢查是否 `n*n` 溢位 ### 線性篩法 > 線性是指演算法時間 $T(n) = O(n)$ 與 Sieve of Eratosthenes 相反:刪所有質數的倍數,就是刪**所有數的質數倍** 而刪質數倍的過程能簡單的判斷是否能及早收手,讓計算效率大幅提升! ```cpp vector<bool> is_p(maxn, true); is_p[1] = false; for (int n = 2; n < maxn; n++) { if (is_p[n]) prime.push_back(n); for (int p: prime) { if (p*n >= maxn) break; // 超出篩檢範圍 is_p[p*n] = false; if (n%p == 0) break; } } ``` #### 存在性 :::info 數到 $n$ 以前,$n$ 就已經被判定過**是否**為質數了。 ::: 根據[唯一分解定理](#質因數分解),$n = p_1 \cdot p_2 \cdot\cdots\cdot p_k$,其中 $p_1\le p_2\le\cdots\le p_k$ 令 $x_1 = p_2 \cdot\cdots\cdot p_k$,由於 $x_1<n$,所以必然有 $p_1$ 與 $x_1$ 相乘得到 $n$ > $p_1$ 就是透過第二層迴圈得來的 > $x_1$ 是在第一層迴圈數到 $n$ 以前曾數到的數字 就算 $p_1$ 整除 $x_1$, `if (n%p == 0) break;` 是之後才執行的,一定能湊得 $n$ #### 唯一性 :::info 任意 $n$ 只會被判定過**恰好一次** ::: 令 $x_i = {n \over p_i}$, $i \in \{1, 2,\cdots, k\}$ $n$ **只**會被數到 $x_1$ 時給**判定**,同樣用 $p_1$ 湊得 $n = p_1\cdot x_1$ 對於 $x_i \not= x_1$, 在數到 $p_1$ 時 $x_i$ 就被 $p_1$ 給整除了,因而不會湊得 $p_i \cdot x_i$ > $p_i$ 根據假設有 $p_i > p_1$ ## 快速冪(Exponentiating by Squaring) 主要目標:給定兩個非負整數 $n, x$,算出 $n^x$。 最直覺的就是跑 $x$ 次迴圈把答案算出來,時間複雜度為 $O(x)$, 但是如果 $x$ 的值非常大(例如 $10^9$),這個方法就很可能會造成程式超時。 > 因為當次方數大的時候非常容易就會超出 `long long` 可以表示的範圍, > 所以通常都會搭配[同餘運算](#同餘運算(Modular-Arithmetic))來使用。 這邊介紹的快速冪可以將時間複雜度降到 $O(\log_2 x)$。 我們可以將 $x$ 分成多個不同的二的冪次的和,例如 $13(1101_2)$ 可以表示為 $1 + 4 + 8$, 則 $n^{13}$ 就可以寫成 $n^{(1 + 4 + 8)} = n \times n^4 \times n^8$,可以看到我們只要去看 $x$ 中哪個位元為 $1$, 將所有這些位元對應的 $n$ 的次方乘起來就是我們要的答案。 下面是範例程式: ```cpp! int pow(int n, int x) { int ans{1}; for (; x; x >>= 1) { if (x & 1) ans *= n; n *= n; } return ans; } ``` > 要注意到因為這邊只是解釋概念, > 因此上面的程式碼只能用來解釋演算法的流程, > 直接拿來執行的話很容易出現整數溢位的錯誤 上面的範例程式中,我們使用 `for` 迴圈來依序從 LSB 檢查到 MSB [^bit_numbering], 並且透過 $n$ 每次都乘以自己來快速得到每個位元對應的 $n$ 的次方。 ### 例題 #### 費式數列 :::info 定義 $f$ 為一個費氏數列,$f_0 = f_1 = 1$, 給定一個整數 $n$,請算出 $f_n$ 與 $10^9 + 7$ 取餘數後的值。 * $0 \leq n \leq 10^{18}$ ::: 這題因為 $n$ 很大,所以不能單純的用動態規劃求解。 我們先將 $f_0, f_1$ 寫成一個 $2 \times 1$ 的矩陣 $F = \left[ \begin{matrix}f_0 \\ f_1 \end{matrix}\right]$,並且定義一個 $2 \times 2$ 的矩陣 $A = \left[ \begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix} \right]$。 可以發現 $AF = \left[\begin{matrix} f_1 \\ f_0 + f_1 \end{matrix}\right] = \left[\begin{matrix} f_1 \\ f_2 \end{matrix}\right]$,而 $A^2F = A(AF) = \left[\begin{matrix} f_2 \\ f_3 \end{matrix}\right]$,也就是 $A^nF = \left[\begin{matrix} f_n \\ f_{n + 1} \end{matrix}\right]$, 所以我們可以先將 $A^n$ 利用快速冪求出來後再與 $F$ 相乘就可以得到 $f_n$。 ```cpp! #include <iostream> #include <array> using namespace std; constexpr int mod{1000000007}; void mtx_mul(array<array<long long, 2>, 2> a, array<array<long long, 2>, 2> b, array<array<long long, 2>, 2> &res) { for (int i{0}; i < 2; ++i) for (int j{0}; j < 2; ++j) { res[i][j] = 0; for (int k{0}; k < 2; ++k) { res[i][j] += a[i][k] * b[k][j] % mod; res[i][j] %= mod; } } } int main(void) { long long n; cin >> n; array<array<long long, 2>, 2> A{0, 1, 1, 1}, ans{1, 0, 0, 1}; // ans 初始值為單位矩陣 for (; n; n >>= 1) { if (n & 1) mtx_mul(ans, A, ans); // ans *= A mtx_mul(A, A, A); // A *= A } cout << (ans[0][0] + ans[0][1]) % mod << '\n'; return 0; } ``` ## 練習題 | Problem | tags | |:-:|:-| | [X-Magic Pair](https://codeforces.com/contest/1612/problem/D) | `gcd` | | [d271: 11582 - Colossal Fibonacci Numbers!](https://zerojudge.tw/ShowProblem?problemid=d271) | `exponentiating by squaring` | | [CODEFORCES 1151C Problem for Nazar](https://codeforces.com/contest/1151/problem/C) | `modulus` | | [CODEFORCES 1165E Two Arrays and Sum of Functions](https://codeforces.com/contest/1165/problem/E) | `modulus`, `sort` | | [ZERO JUDGE a289 Modular Multiplicative Inverse](https://zerojudge.tw/ShowProblem?problemid=a289) | `modulus` | | [NCKU OJ 22 爬樓梯 k](https://oj.leo900807.tw/problems/22) | `dp`, `modulus`, `combinatorics` | | [CODEFORCES 1033B Square Difference](https://codeforces.com/problemset/problem/1033/B) | `prime` | | [ZERO JUDGE a007 判斷質數](https://zerojudge.tw/ShowProblem?problemid=a007) | `prime`, `primality test` | | [CODEFORCES 1165D Almost All Divisors](https://codeforces.com/contest/1165/problem/D) | `number theory` | | [CODEFORCES 1114C Trailing Loves (or L'oeufs?)](https://codeforces.com/contest/1114/problem/C) | `brute force`, `number theory` | | [LeetCode 952 Largest Component Size by Common Factor](https://leetcode.com/problems/largest-component-size-by-common-factor/) | `dsu`, `math` | | [UVa OJ 543 Goldbach’s Conjecture](https://uva.onlinejudge.org/external/5/543.pdf) | `prime`, `sieve` | | [UVa OJ 10140 Prime Distance](https://uva.onlinejudge.org/external/101/10140.pdf) | `prime`, `sieve` | ## Reference * [《提升程式設計的解題思考力 - 國際演算法程式設計競賽訓練指南》Chapter 02 數學基礎](https://hackmd.io/@XDEv11/S1sSbRu3d) * [貝祖定理](https://zh.wikipedia.org/wiki/%E8%B2%9D%E7%A5%96%E7%AD%89%E5%BC%8F) [^bit_numbering]: [Bit significance and indexing](https://en.wikipedia.org/wiki/Bit_numbering#Bit_significance_and_indexing) [^coprime]: 最大公因數為 $1$ [^prod]: $\prod$ 是**連乘**符號,也就是說 $\prod\limits_{k=1}^n m_k = m_1\times m_2\times \cdots \times m_n$ [^probability_a]: 此機率是根據測試的次數計算,每次測試會有不同的 $a$