owned this note
owned this note
Published
Linked with GitHub
---
type: slide
slideOptions:
transition: fade;
---
<style>
.reveal .slides {
text-align: left;
font-size:30px;
}
</style>
# Number Theory
Introduction to Competitive Programming
2024/4/17
藍白
---
* 質因數分解
* 鴿籠原理 :dove_of_peace:
* 模運算 & 模逆元
* 歐幾里得算法(輾轉)
* 拓展歐幾里得(求逆元，求方程解)
* 中國剩餘定理(求方程解)
* 排容定理(求聯集)
---
Disclaimer: 任何在這講義上的程式碼，請先用一堆題目驗證正確性，方可在比賽中使用
Jakao 推薦用[這網站](https://judge.yosupo.jp/)驗證
---
## 質數檢驗 & 質因數分解
~~上過了，所以跳過。~~
並沒有，我們來上點好玩的
以下東西我們只負責用，不負責證明，證明請工五松苑左轉鴻經館或看 OI wiki
----
## 質數檢驗大絕：Miller–Rabin
時間複雜度 $O(k \cdot lg^3(n))$，k 為 k 論測試，n 為待測數字
~~利用快速傅里葉轉換等技術可以將時間複雜度壓一個次方，但我不想研究~~
----
## Templete
* 以下程式碼有過 Zerojudge a007（~~指標性的驗質數題目~~）
* 可以驗證出 $2^{64}$ 以內的所有質數
```cpp=
LL modMul(LL a, LL b, LL m){
return __int128{a} * b % m;
}
LL pow(LL a, LL b, LL m){LL ret = 1;
for (; b; a = modMul(a, a, m), b >>= 1)
if (b % 2) ret = modMul(ret, a, m);
return ret;
}
bool isPrime(LL n){
LL sprp[7] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
if(n == 1 || (n & 1) == 0) return n == 2;
int u = n - 1, t = 0; for(; u % 2 == 0; t++) u >>= 1;
for(int i = 0; i < 7; i++){ LL a = sprp[i] % n;
if(a == 0 || a == 1 || a == n - 1) continue;
LL x = pow(a, u, n); if(x == 1 || x == n-1) continue;
for(int j = 1; j < t; j++){ x = modMul(x, x, n);
if (x == 1) return 0; if (x == n - 1) break;
}
if(x == n - 1)continue; return 0;
}
return 1;
}
```
----
## 質因數分解大絕：Pollard's Rho
* 目前暫時沒有用過
```cpp=
LL pollard_rho(LL n){
//pre-define f = (x * x + 1) % mod
if(!(n&1)) return 2;
while(1){
LL y = 2, x = rand()%(n-1) + 1, res = 1;
for(int sz = 2; res == 1; sz *= 2){
for(int i = 0; i < sz && res <= 1; i++){
x = f(x, n);
res = __gcd(abs(x - y), n);
}
y = x;
}
if(res != 0 && res != n) return res;
}
}
```
----
痛苦題目：2022 TOPC pE: Etched Emerald Orbs
題意：
給你 $3 < k < 4 \cdot 10^{18}$，輸出 x, y 使得 $2/k = 1/x + 1/y$
其中 $1 < x < y < 2^{125}$，且要求有多組解時，$x+y$ 總和最小
----
* 官方參考解答
* 純欣賞用
```cpp
#include <bits/stdc++.h>
using namespace std;
#define O ostream
O& operator << (O &out, __int128_t v) {
O::sentry s(out);
if (s) { __uint128_t uv = v < 0 ? -v : v;
char buf[128], *d = end(buf);
do { *(--d) = uv % 10 + '0'; uv /= 10;
} while (uv != 0);
if (uv < 0) *(--d) = '-'; int len = end(buf) - d;
if (out.rdbuf()->sputn(d, len) != len)
out.setstate(ios_base::badbit);
}
return out;
}
#define I istream
I& operator >> (I &in, __int128_t &v) {
string s; in >> s; v = 0;
for (int i = 0 ; i < (int)s.size() ; i++)
if ('0' <= s[i] && s[i] <= '9')
v = 10 * v + s[i] - '0';
return in;
}
using LL = long long;
LL modMul(LL a, LL b, LL m) {
return __int128{a} * b % m;
}
LL pow(LL a, LL b, LL m) { LL ret = 1;
for (; b; a = modMul(a, a, m), b >>= 1)
if (b % 2) ret = modMul(ret, a, m);
return ret;
}
bool isPrime(LL n) {
LL sprp[7] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
if (n == 2) return true; if (n < 2 || n % 2 == 0) return false;
LL t = 0, u = n - 1; for (; u % 2 == 0; t++) u >>= 1;
for (int i = 0 ; i < 7 ; i++){ LL a = sprp[i] % n;
if (a == 0 || a == 1 || a == n-1) continue;
LL x = pow(a, u, n); if (x == 1 || x == n-1) continue;
for (int j = 1 ; j < t ; j++){ x = modMul(x, x, n);
if (x == 1) return false; if (x == n - 1) break;
}
if (x == n - 1) continue; return false;
}
return true;
}
inline LL f(LL x , LL mod) {
return (modMul(x, x, mod) + 1) % mod;
}
inline LL pollard_rho(LL n) {
if(!(n&1)) return 2;
while(true) {
LL y = 2 , x = rand() % (n - 1) + 1 , res = 1;
for(int sz = 2; res == 1; sz *= 2) {
for(int i = 0; i < sz && res <= 1; i++) {
x = f(x , n);
res = __gcd(abs(x - y) , n);
}
y = x;
}
if (res != 0 && res != n) return res;
}
}
void factorize(LL n, auto &ans) {
if (isPrime(n))
ans.push_back(n);
else {
auto p = pollard_rho(n);
factorize(p, ans);
factorize(n / p, ans);
}
}
// 枚舉因數
// int i 指第幾個質因數，dfs 內部的兩個 dfs 分別代表該質因數的 0 次方與 j 次方的情形
void dfs(int i, auto &pool, __int128 x, auto &factors) {
if (i == pool.size()) {
factors.push_back(x);
return ;
}
dfs(i + 1, pool, x, factors);
for (int j = 0; j < pool[i].second; j++) {
x *= pool[i].first;
dfs(i + 1, pool, x, factors);
}
}
int main() {
LL k; cin >> k;
assert(3 <= k and k <= 4'000'000'000'000'000'000);
auto p = vector<LL>{};
factorize(k, p);
map<LL, int> cnt;
for (auto &v : p) {
cnt[v] += 2;
}
// 將 map<LL, int> 單純轉成 vector<pair<LL, int>>
auto pool = vector<pair<LL, int>>(cnt.begin(), cnt.end());
auto factors = vector<__int128>{};
dfs(0, pool, 1, factors);
// 因數由大到小排序
sort(factors.rbegin(), factors.rend());
auto ans_x = __int128{k} * k;
auto ans_y = __int128{k} * k;
for (auto &p : factors) {
if ((p + k) % 2 != 0)
continue;
auto x = (p + k) / 2;
auto y = __int128{k} * k / p + k;
if (y % 2 != 0)
continue;
y /= 2;
if (x >= y)
continue;
if (x + y < ans_x + ans_y) {
ans_x = x;
ans_y = y;
}
}
cout << ans_x << ' ' << ans_y << '\n';
}
```
---
## 鴿籠原理
----
亦稱抽屜原理。
常用於證明存在性證明(有兩個或以上)和求最壞情況(最差會有一個)下的解
經典問題就是鴿籠問題，抽象一點地說，將 $n+1$ 個物體劃分為 $n$ 個組別，則至少會有一個組別有兩個或以上的物體。
若今天推廣到 $k$ 個物體，分成 $n$ 組，則至少會有一組有 $\lceil \frac k n \rceil$個(或以上)個物體。
----
## [實際應用在題目](https://vjudge.net/problem/UVA-11237)
關鍵點為製造抽屜(製造籠子)
ex : 給出 $n$ 個數字，和整數 $m(1 \leq m \leq n-1)$，問說是否可以選出至少兩個數字都同餘 $m$ 。
解題思路 : $n > m$，製造 $m$ 個籠子
----
## [Zerojudge c268](https://zerojudge.tw/ShowProblem?problemid=c268)
題意：給你 n 段長度，請問這些長度能不能構成一個三角形？
其實這題不大鴿籠，但要表達的是那種「塞不下」的感覺
利用類似 Greedy 的手段可以知道，只要超過 $n > 45$ 就必有一個三角形存在，$n \le 45$ 就排序之後一個迴圈跑 $O(n)$ 個 a[i], a[i+1], a[i+2] 就好
---
## 模運算 & 模逆元
----
## 模運算
整數在模空間底下計算，要讓計算結果在 $[0, m-1]$ 之間
加減乘可以直接做，但除法會發生甚麼事情呢?
(此為常用技巧)
$(a + b) \mod m = (a+b)\ \%\ m$
$(a - b) \mod m = ((a-b)\ \%\ m + m)\ \%\ m$
$(a \cdot b) \mod m = (a \cdot b)\ \%\ m$
$(a\ /\ b) \mod m =\ ?$
----
## [模逆元](https://atcoder.jp/contests/abc154/tasks/abc154_f)
發現除法會錯(也避免負數)之後，於是我們知道，需要 模 "逆元" !
在模 $m$ 意義下，對於一個整數 $u$，如果存在一個整數 $v$，使得 $u \cdot v≡1(mod\ m)$，則 $v$是 $u$ 的逆元，記作 $u^{−1}$
因此對於 $(a\ /\ b)\mod m$
如果找到 $b$ 的模逆元就可以把除法運算改成乘法了，也就是尋找**正整數** $v = b^{-1}$ 使得
$(a\ /\ b)\mod m = (a \cdot b^{-1})\mod m = (a \cdot v) \mod m$
----
## 逆元性質以及證明
* 存在唯一性
模p意義下，對於a而言，如果其存在逆元，則逆元只有一個$a′ = a^{−1}$。
可使用反證法或演譯法證明。
* 完全積性函數
也就是 $(a \cdot b)^{−1}≡a^{−1} \cdot b^{−1}(mod \ p)$
考慮有 $a \cdot a^{−1}≡1(mod \ p)$，和$b \cdot b^{−1}≡1(mod \ p)$，乘到一起就是$(a \cdot b)(a^{−1}×b^{−1})≡1(mod \ p)$，所以$(a \cdot b)^{−1}≡a^{−1}×b^{−1}(mod \ p)$
[證明參考](https://www.cnblogs.com/ac-evil/p/12680705.html)
於是就可以發現除法可以變成$(a\ /\ b) \equiv (a\ \times\ b^{-1}) \mod m$
----
## 費馬小定理
如果 $a, m$ 互質的情況下
$a^m \equiv a (\mod m)$
$\to$
$a^{m-1} \equiv 1 (\mod m)$
$\to$
$a^{m-2} \equiv a^{-1} (\mod m)$
因此求出 $a^{m-2}$ 即為 $a$ 的模逆元
**而當 a, m 不互質 ( gcd $\ne$ 1) 的情況下模逆元不存在**
----
因此在 $b$ 與 $m$ 互質的情況下
要求出 $a / b$ 可以求出
$(a\ /\ b) \equiv (a\ \times\ b^{m-2}) \mod m$
而 $b^{m-2}$ 可以透過快速冪求出
通常需要取模的題目給的 $m$ 都會是很大的質數 ($10^9+7, 998244353$) 等
所以通常 $b$ 會跟 $m$ 互質
---
# 下課 :bread:
---
## 歐幾里得算法 (輾轉相除法)
----
計算兩個數字 $a$ , $b$ 的最大公因數 $gcd(a,b)$
同時在程式實作裡，我們會使用輾轉相除法進行實作。
即 $gcd(a,b)=gcd(b,a \ mod \ b)$ $a$ , $b$ 非負且 $a\ge b$
輾轉相除法求解最大公因數並不是唯一方法，但相當快速，時間複雜度為 $O(lg(max(a, b)))$
----
## 具體實作
```cpp
int gcd(int a, int b) {
int r;
while (b) {
r = a%b, a = b, b = r;
}
return a;
}
```
那其實也可以寫成
```cpp
int gcd(int a,int b) {
return b == 0 ? a : gcd(b,a % b);
}
```
----
當然也可以直接呼叫預設函式
since C++14
```cpp
#include<algorithm>
__gcd(a,b);
```
since C++17
```cpp
#include<numeric>
gcd(a,b);
```
---
## 拓展歐幾里得--EXGCD
----
## 裴蜀定理
對於不定方程$ax+by=m$的整數解，其有解的充要條件為 $gcd(a,b)|m$
$a, b \in N, \; gcd(a,b) = d$
$\forall x\forall y \; d\:|\:a x + b y$
$\exists x\exists y \; a x + b y = d$
----
因為輾轉相除法除到最後餘數為 0，在這，我們設 $r_{k+1}=0$，那麼， $r_k$ 就是 $a$ 和 $b$ 的最大公因數，即 $r_k = d$。將 $r_k = d$帶入上式中可以移項得到
$d=m_1 \cdot r_{k−2}+n_1 \cdot r_{k−3}$
----
由於其中會用到的數字都是整數，在推論到最後會發現式子會變成
$d=m_k \cdot a+n_k \cdot b$
於是得證，$ax + by = d$ 有整數解
最後，需要證明方程$ax+by=z$，只有滿足 $gcd (a,b)∣z$ ，方程才有整數解。
----
## 裴蜀定理在歐幾里得上的應用
於是，知道$ax+by=m$的整數解其有解的充要條件為 $gcd(a,b)|m$ 後，就可以回到正題，要如何在已知 $a$ , $b$ 的情況求解式子裡的$x$ , $y$
設當前的式子為 $ax_1 + by_1 = gcd(a,b)$，根據歐幾里得算法，存在式子
$bx_2+(a \ mod \ b)y_2 = gcd (b,a mod b)$
$∵ a \ mod \ b = a − \lfloor \frac a b \rfloor \cdot b$
$∴ ax_1 + by_1 = bx_2 + (a− \lfloor \frac a b \rfloor \cdot b)y_2$
$ax_1 + by_1 = bx_2 + ay_2 − \lfloor \frac a b \rfloor \cdot b \cdot y_2$
$ax_1 + by_1 = ay_2 + b(x_2 − \lfloor \frac a b \rfloor \cdot y_2 )$
$∴ x_1 = y_2 , y_1 = x_2 − \lfloor \frac a b \rfloor \cdot y_2$
```cpp
int exgcd(int a,int b,long long &x,long long &y) {
if(b == 0){x=1,y=0;return a;}
int now=exgcd(b,a%b,y,x);
y-=a/b*x;
return now;
}
```
----
exgcd 還有甚麼用途呢?可以用來求逆元。
求 $a$ 在模 $p$ 意義下的逆元
倘若使用費馬小定理則需要保證p為質數，使用 exgcd 則不需要。
設 $x$ 為 $a$ 在模 $p$ 意義下的逆元，那麼滿足式子：
$ax \equiv 1 (mod \ p)$
那麼有：
$ax + my = 1$
然后用 exgcd 求出 $x$ 即可
```cpp
long long inv(long long a,long long m){
long long x,y;
long long d=exgcd(a,m,x,y);
if(d==1) return (x+m)%m;
else return -1; //-1為無解
}
```
----
## Problem
### NCPC 2020 final
----
```
input
3
2 3 5
4 4 3
13 5 7
output
2
1
5
```
----
只需要移項完式子後使用三分搜搜索 $x$ 及 $y$ 即可，會發現他的 $x$ 值以及 $y$ 值的圖形是為一個斜線(直線)，所以只需要使用三分搜查詢到那兩個點，就可以確定答案。
---
## 中國剩餘定理 (Chinese Remainder Theorem, CRT)
:bread:
----
也是韓信點兵那些的經典例題，應該要對這些東西很熟悉
中國剩餘定理:設正整數 $m_1$ , $m_2$ …… $m_k$ 兩兩互質，則同餘方程組求解 $x$
$x \equiv a_1 \pmod {m_1}$
$x \equiv a_2 \pmod {m_2}$
$x \equiv a_3 \pmod {m_3}$
.
.
.
$x \equiv a_{k-2} \pmod {m_{k-2}}$
$x \equiv a_{k-1} \pmod {m_{k-1}}$
$x \equiv a_{k} \pmod {m_{k}}$
----
那如何求 $x$ 呢?
過程:
1. 計算所有模數的積 $n$
2. 對於第 $i$ 個方程：
計算 $M_i=\frac{n}{m_i}$
計算 $M_i$ 在模 $m_i$ 意義下的 逆元 $M_i^{-1}$
方程組在模 $n$ 意義下的唯一解為：$x=\sum_{i=1}^k a_i \cdot M_i \cdot M_i^{-1} \pmod n$
----
## 實作方式 會使用到模逆元，也就是會用到拓展歐幾里得
```cpp
#define LL long long
LL exgcd(LL a,LL b,LL &x,LL &y){
if(!b){
x = 1;
y = 0;
return a;
}
int now=exgcd(b, a % b, y, x);
y -= a / b * x;
return now;
}
LL CRT(LL k, LL* a, LL* r) {
LL n = 1, ans = 0;
for (LL i = 1; i <= k; i++) {
n = n * r[i];
}
for (LL i = 1; i <= k; i++) {
LL m = n / r[i], b, y;
exgcd(m, r[i], b, y);
ans = (ans + a[i] * m * b % n) % n;
}
return (ans % n + n) % n;
}
```
----
## Problem
### NCPC 2018 final
----
```
input
3
2 4 1
4 5 9
1 2 4
2 5 7
0 0 1
2 5 127
output
154
67
890
```
----
跟範例程式碼差不多，只需要稍微確認輸入的順序，以及好好的初始化，然後套入程式碼即可。
---
## 排容定理 (容斥原理)
如果我們使用$A , B , C$ 去代表三個集合，則元素總個數會是 $|A \cup B \cup C|$
那如果我們要計算這個，則需要扣除重複的部分，也就是$|A \cap B| , |B \cap C| , |C \cap A|$
那會發現我們又不小心多刪除了，於是需要加回$|A \cap B \cap C|$
也就是
$|A \cup B \cup C| = |A| + |B| + |C| - |A \cap B| - |A \cap C| - |B \cap C| + |A \cap B \cap C|$
----
在 $N$ 個集合 $A_1, A_2,..., A_n$裡，如果要算出這 $N$ 個集合的聯集
$|\bigcup\limits_{i=1}^{N} A_i| =$
$\sum\limits_{i=1}^{N}|A_i| - \sum_{1 \le i < j \le N}|A_i\cap A_j| + \sum_{1 \le i < j < k \le N}|A_i\cap A_j\cap A_k| - ...$
$+ (-1)^{N-1}|A_1\cap ...\cap A_N|$
一個集合 減去 兩個集合的交集 加上 三個集合的交集 ... 直到算出所有可能的交集為止。
----
使用之前教過的，大家也看過的圖片當作例子
----
## 二項式定理
要如何知道每個集合之間的聯集要加幾次，上圖片
會發現可以充分利用帕斯卡三角形的性質，對於每一個項數去進行補充或削減，直到每個元素的出現次數皆為$1$時即可。
----
## 例題 [CSES 2185. Prime Multiples](https://cses.fi/problemset/task/2185)
給你數字 $N, M$
給 $N$ 個相異質數
問在 $1\sim M$ 之間有幾個數字為至少一個給定的質數的倍數
- $1 \le N \le 20$
- $1 \le M \le 10^{18}$
### sample
```
N = 3, M = 15
質數為 {2, 5, 7}
10
(2,4,5,6,7,8,10,12,14,15)
```
----
因此 $N$ 個集合為 第 $i$ 個質數的倍數
要求出 $1\sim M$ 之間的數字所有集合的聯集
以範例來說
$|2|$ + $|5|$ + $|7|$ - $|2\cap 5|$ - $|2\cap 7|$ - $|5\cap 7|$ + $|2\cap 5\cap 7|$
因此答案為 7 + 3 + 2 - 1 - 1 - 0 + 0 = 10
----
範例程式碼
```cpp
for (int i = 1; i < 1<<k; i++) {
int x = -1;
if (__builtin_popcount(i)&1) x = 1;
int y = n;
int z = 1;
for (int j = 0; j < k; j++) {
if (i>>j&1) {
if (z >= n/a[j] + 1) {
z = INF;break;
}
z = z*a[j];
}
}
y /= z,ans += x*y;
}
```
---
### ICPC 2019 Taipei-Hsinchu regional pM: DivModulo
題意：
若 $n= m \cdot d^x$，且 d, m 互質，定義 $n\;dmod\;d = m\;mod\;d$
給定 $M, N, D$，求 $C(M,N)\;dmod\;D$，其中$N \le M\le 4 \cdot 10^{18}$, $2 \le D \le 1.6 \cdot 10^7$
**因為我也不會，所以這題純觀賞用**
---
### 如果你實在很愛數論
https://oi-wiki.org/math/
Google Drive
Gist
Clipboard
Sign in with Wallet
