:+1: [2020 競技程設教材 HackMD Book](https://hackmd.io/@nckuacm/ryLIV6BYI) :+1:

2019 Week 15: Number theory & Calculation
=
>數大便是美

# 模反元素

數字們同餘 $n$ 後，對於加法、減法、乘法，其性質都與以往差不多，但對數字做**除法**卻不盡人意
> 注意，在同餘運算中只會有**整數**，有理數無理數等其他數不會出現

**反元素**指的是元素與其**運算**後為**單位元素**的元素
例如：
- 整數 $a$ 的**加**法反元素為 $-a$
- 整數 $a$ 的**乘**法反元素為 $a^{-1}$

而元素 $a$ 的**模 $n$ 反元素**為 $a^{-1}$ 滿足
$$ a \cdot a^{-1} \equiv 1 \mod n$$

根據 [Bezout's Thm](https://hackmd.io/4nPaMD6wR02JqsBCEgNimw?view#Bezout%E2%80%99s-Theorem)
- $ax + ny = 1$ 在模 $n$ 後有 $ax \equiv 1 \mod n$
- $ax + ny \not= 1$ 在模 $n$ 後有 $ax \not\equiv 1 \mod n$

也就是說 $\gcd(a, n) = 1$ (互質)，$a$ 才擁有模 $n$ 反元素

## 求出模反元素
- 當 $n$ 為質數時：

根據[費馬小定理](https://hackmd.io/4nPaMD6wR02JqsBCEgNimw?view#%E8%B2%BB%E9%A6%AC%E5%B0%8F%E5%AE%9A%E7%90%86)有
$$a \cdot a^{n-2} \equiv 1 \mod n$$
所以 $a^{n-2}$ 是 $a$ 的模反元素，可用[快速冪](#快速冪) 在 $O(\log n)$ 算出

- 當 $n$ 非質數時：

根據 Bezout's Thm 有
$$
ax + ny = 1 \Rightarrow a \cdot x \equiv 1 \mod n
$$
可用[擴展歐幾里得演算法](https://hackmd.io/4nPaMD6wR02JqsBCEgNimw?view#Extended-Euclidean-algorithm)找到 $x = a^{-1}$

#### 練習：
[ZERO JUDGE a289 Modular Multiplicative Inverse](https://zerojudge.tw/ShowProblem?problemid=a289)

# 質數判斷
若數 $n > 1$ 能被 $1$ 與 $n$ 以外的數整除，則 $n$ 為合數

bool is_p = true; // is_p := is it prime?
for (int i = 2; i <= sqrt(n); i++)
  if (n%i == 0) is_p = false;

#### 練習：
[CODEFORCES 1033B Square Difference](https://codeforces.com/problemset/problem/1033/B)

# 質因數分解
根據唯一分解定理，任何合數都能被分解成一些質數的積
> 算術基本定理

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});

#### 練習：
[CODEFORCES 1165D Almost All Divisors](https://codeforces.com/contest/1165/problem/D)
[CODEFORCES 1114C Trailing Loves (or L'oeufs?)](https://codeforces.com/contest/1114/problem/C)
[LeetCode 952 Largest Component Size by Common Factor](https://leetcode.com/problems/largest-component-size-by-common-factor/)
# 質數篩檢
>數質數可以有效安撫緊張的情緒

## Sieve of Eratosthenes
其精神是將質數的**倍數**都設為**非質數**


fill(is_p, 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` 溢位

#### 練習：
[UVa OJ 543 Goldbach’s Conjecture](https://uva.onlinejudge.org/external/5/543.pdf)
[UVa OJ 10140 Prime Distance](https://uva.onlinejudge.org/external/101/10140.pdf)
\* [ZERO JUDGE a007 判斷質數](https://zerojudge.tw/ShowProblem?problemid=a007) [^萬惡a007]

## 線性篩法

fill(is_p, 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 = p_1 \cdot p_2 \cdot ...\cdot p_k$，其中 $p_1\le p_2\le...\le p_k$

令 $x_1 = p_2 \cdot ...\cdot p_k$，由於 $x_1<n$，所以必然有 $p_1$ 與 $x_1$ 相乘得到 $n$
> $p_1$ 就是透過第二層迴圈得來的

就算 $p_1$ 整除 $x_1$， `if (n%p == 0) break;` 是之後才執行的，一定能湊得 $n$

### 唯一性
:::info
任意 $n$ 只會被判定過**恰好一次**
:::

令 $x_i = {n \over p_i}$， $i \in \{1, 2,..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$

# 算術
對於**大數字**的運算，普通的做法不夠快，因此接下來將介紹快速的**乘法**及**冪**運算

## Gauss's complex multiplication algorithm
對於**複數** $a+bi$ 與 $c+di$ **相乘**
$$
(a+bi) \cdot (c+di) = (ac-bd) + (bc+ad)i
$$
相較於 $4$ 次乘法 $ac, bd, bc, ad$

可只用 $3$ 次乘法 $\begin{cases}k_1 = c(a+b)\\ k_2=b(c+b)\\ k_3=a(d-c)\end{cases}$ 完成複數乘法：

$$
\text{where } \begin{cases}
\begin{split}
(ac - bd) &= ac + \mathbf{bc} - \mathbf{bc} - bd &= \underline{c\cdot (a+b)} - \underline{b\cdot (c+b)} =k_1-k_2\\
(bc + ad) &= \mathbf{ac} + bc + ad -\mathbf{ac} &= \underline{c\cdot (a+b)} + \underline{a\cdot (d-c)} =k_1+k_3
\end{split}
\end{cases}
$$
因此複數相乘 $$(a+bi) + (c+di) = (k_1-k_2) + (k_1+k_3)i$$

## Karatsuba algorithm
對於**整數** $\begin{cases}\begin{split}x &= am+b\\y &= cm+d \end{split}\end{cases}$ **相乘**
> 受到複數乘法啟發，將整數 $x$ 分割成兩數 $a, b$

$$
x\cdot y = ac\cdot m^2 + (ad+bc)\cdot m + bd
$$

能用 $3$ 次乘法 $\begin{cases}z_1=ac\\ z_2=bd\\ z_3=(a+b)(c+d)\end{cases}$ 就完成運算：
$$
\begin{split}
(ad+bc) &= \mathbf{ac}+ad+bc+\mathbf{bd}-\mathbf{ac}-\mathbf{bd} \\
&=\underline{(a+b)(c+d)}-\underline{ac}-\underline{bd} \\
&= z_3-z_1-z_2
\end{split}
$$

因此整數相乘 $$x\cdot y=z_1m^2 + (z_3-z_1-z_2)m+z_2$$

#### 範例 $\begin{cases}\begin{split}12345 &= 12 &\cdot 1000 &+345\\6789 &= 6 &\cdot 1000 &+789 \end{split}\end{cases}$相乘：

首先 $\begin{cases}z_1=12\cdot 6 &= 72\\ z_2=345\cdot 789 &= 272205\\ z_3=(12+345)(6+789) = 357\cdot 795&=283815\end{cases}$

則 $\begin{split}12345 \cdot 6789 &= 72 \cdot 1000^2 + (283815-72-272205)\cdot 1000+272205 \\
&=72000000+11538000+272205 \\ &= 83810205
\end{split}$

> 似乎沒有比較快欸
> 那是因為你**想**的乘法，仍然是 $O(n^2)$ 的[直式乘法](https://en.wikipedia.org/wiki/Multiplication_algorithm#Long_multiplication)

### 分治法
對於上述演算法，凡是遇到乘法運算，都使用同樣的演算法
並且對於數字的分割，總是分割成均等的**兩半**
> 例如上面 $6789$ 應分成 $67$ 和 $89$
> $6789 = 67\cdot 100 + 89$

int k(int x, int y) {
  if(x < 10 || y < 10) return x*y;

  int len = min(log10(x), log10(y));
  int m = pow(10, len/2 + 1);

  auto [a, b] = div(x, m); // since c++17
  auto [c, d] = div(y, m);

  int z1 = k(a, c);
  int z2 = k(b, d);
  int z3 = k(a+b, c+d);

  return z1*m*m + (z3-z1-z2)*m + z2;
}

時間成本 $T(n)$ **約**為
$$
\begin{split}
T(n) &= 3 \cdot T(\lceil n/2 \rceil) + c_n\\
T(1) &= 1 + c_1
\end{split}
$$
此處 $n$ 為整數長度，$c_k$ 為**加法**運算的成本
複雜度為 $O(3^{\log_2n}) = O(n^{\log_23})$

## 快速冪
求 $a$ 的 $n$ 次**冪** $x =a^n$：

int x = 1;
while (n--) x *= a;

基於 D&C，若 $n$ 為**偶數**則 $a^n = a^{n\over 2} \times a^{n\over 2}$

- top-down

int fast_exp(int a, int n) {
  if (n == 1) return a;
  
  int x = fast_exp(a, n/2);
  return (x*x) * (n&1? a : 1); // 檢查是否為奇數
}

並且 $a^{n+m} = a^n \times a^m$
- bottom-up

int x = 1;
while (n) {
  if (n&1) x *= a;
  a *= a;
  n >>= 1;
}

複雜度從樸素的 $O(n)$ 改進至 $O(\log_2n)$

#### 練習：
[UVa OJ 374 Big Mod](https://uva.onlinejudge.org/external/3/p374.pdf)
[UVa OJ 10006 Carmichael Numbers](https://uva.onlinejudge.org/external/100/10006.pdf)
[CODEFORCES 615D Multipliers](https://codeforces.com/contest/615/problem/D)

## 矩陣快速冪
對於**方陣**乘法，也能沿用快速冪的想法
例如求第 $n$ 項的費式數 $f_n$：

$$
\begin{bmatrix}
1 &1 \\
1 &0 \\
\end{bmatrix}^n
\cdot
\begin{bmatrix}
1 \\
0
\end{bmatrix} =
\begin{bmatrix}
f_{n+1} \\
f_n
\end{bmatrix}
$$

int M[][2] = {{1, 1}, {1, 0}};
int f[] = {1, 0};

while (n) {
  if (n&1) {
    int t[] = {f[0], f[1]};
    f[0] = M[0][0] * t[0] + M[0][1] * t[1];
    f[1] = M[1][0] * t[0] + M[1][1] * t[1];
  }

  int p[][2] = {{M[0][0], M[0][1]}, {M[1][0], M[1][1]}};
  M[0][0] = p[0][0] * p[0][0] + p[0][1] * p[1][0];
  M[0][1] = p[0][0] * p[0][1] + p[0][1] * p[1][1];
  M[1][0] = p[1][0] * p[0][0] + p[1][1] * p[1][0];
  M[1][1] = p[1][0] * p[0][1] + p[1][1] * p[1][1];

  n >>= 1;
}

其複雜度 $O(\log_2 n)$

#### 練習：
[ZERO JUDGE b525 先別管這個了，你聽過turtlebee嗎？](https://zerojudge.tw/ShowProblem?problemid=b525)

[^萬惡a007]: 這根本不是基礎題= =