---
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
其精神是將質數的**倍數**都設為**非質數**(合數)
```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$
Sign in with Wallet
Connect another wallet