2020 競技程設教材 HackMD Book

2020 Week 15: Number theory & Calculation

數大便是美

本週繼續由第十四週主題延伸，
也就是說以下內容或多或少需要最大公因數以及模運算的相關定理。

大數乘法

大數字的乘法，普通的直式乘法

$O(N^2)$ 不夠快，因此將介紹快速的乘法運算

這裡

$N$ 為數字的位數

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{array}{l}{k}_1=c(a+b)\\ k_2=b(c+b)\\ k_3=a(d-c)\end{array}$ 完成複數乘法：

$$\text{where }\{\begin{array}{l}\begin{array}{rlr}(ac-bd)& =ac+\mathbf{bc}-\mathbf{bc}-bd& =\underset{―}{c\cdot (a+b)}-\underset{―}{b\cdot (c+b)}=k_1-k_2\\ (bc+ad)& =\mathbf{ac}+bc+ad-\mathbf{ac}& =\underset{―}{c\cdot (a+b)}+\underset{―}{a\cdot (d-c)}=k_1+k_3\end{array}\end{array}$$
因此複數相乘

$$\begin{array}{rl}(a+bi)\cdot (c+di)& =(ac-bd)+(bc+ad)i\\ & =(k_1-k_2)+(k_1+k_3)i\end{array}$$

乘法操作變少，但加減法操作變多了(？)

Karatsuba algorithm

對於整數

$\{\begin{array}{l}\begin{array}{rl}x& =am+b\\ y& =cm+d\end{array}\end{array}$ 相乘

受到複數乘法啟發，將整數

$x$ 分割成兩數

$a,b$

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

能用

$3$ 次乘法

$\{\begin{array}{l}{z}_1=ac\\ z_2=bd\\ z_3=(a+b)(c+d)\end{array}$ 就完成運算：

$$\begin{array}{rl}(ad+bc)& =\mathbf{ac}+ad+bc+\mathbf{bd}-\mathbf{ac}-\mathbf{bd}\\ & =\underset{―}{(a+b)(c+d)}-\underset{―}{ac}-\underset{―}{bd}\\ & =z_3-z_1-z_2\end{array}$$

因此整數相乘

$$x\cdot y=z_1{m}^2+(z_3-z_1-z_2)m+z_2$$

範例

$\{\begin{array}{l}\begin{array}{rlrl}12345& =12& \cdot 1000& +345\\ 6789& =6& \cdot 1000& +789\end{array}\end{array}$ 相乘：

首先

$\{\begin{array}{ll}{z}_1=12\cdot 6& =72\\ z_2=345\cdot 789& =272205\\ z_3=(12+345)(6+789)=357\cdot 795& =283815\end{array}$

則

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

似乎沒有比較快欸
那是因為你想的乘法，仍然是

$O(n^2)$ 的直式乘法

對於上述演算法，凡是遇到乘法運算(求

$z_1,z_2,z_3$)，都使用同樣的演算法
並且對於數字的分割，總是分割成均等的兩半

例如範例中

$6789$ 應分成

$67$ 和

$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{array}{rl}T(n)& =3\cdot T(\lceil n/2\rceil )+c_n\\ T(1)& =1+c_1\end{array}$$
此處

$n$ 為整數位數，

$c_n$ 為加法運算的成本
複雜度為

$O(3^{{\log }_{2}n})=O(n^{{\log }_{2}3})$

---

求大數冪

樸素的求整數

$a$ 的

$n$ 次冪

$x=a^n$，可以採用連乘：

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

整數平方求冪 (Exponentiating by squaring)

又稱快速冪演算法

基於分治法，

• 若
  $n$ 為偶數則
  $a^n=a^{\frac{n}{2}}\times a^{\frac{n}{2}}$
• 若
  $n$ 為奇數則
  $a^n=a^{\lfloor \frac{n}{2}\rfloor }\times a^{\lfloor \frac{n}{2}\rfloor }\times a$

Top-down

int power(int a, int n) {
  if (n == 1) return a;
  
  int x = power(a, n/2);
  return n&1? x*x*a : x*x;
}

Bottom-up

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

複雜度從樸素的

$O(n)$ 改進至

$O({\log }_{2}n)$

練習：

UVa OJ 374 Big Mod
CODEFORCES 615D Multipliers

矩陣平方求冪

對於方陣乘法，也能沿用快速冪的做法
例如求第

$n$ 項的費式數

$f_n$：

$$\begin{array}{c}{\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^n\cdot \left[\begin{array}{c}1\\ 0\end{array}\right]=\left[\begin{array}{c}{f}_{n+1}\\ f_n\end{array}\right]\end{array}$$

Bottom-up

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 t[][2] = {{M[0][0], M[0][1]}, {M[1][0], M[1][1]}};
  M[0][0] = t[0][0] * t[0][0] + t[0][1] * t[1][0];
  M[0][1] = t[0][0] * t[0][1] + t[0][1] * t[1][1];
  M[1][0] = t[1][0] * t[0][0] + t[1][1] * t[1][0];
  M[1][1] = t[1][0] * t[0][1] + t[1][1] * t[1][1];

  n >>= 1;
}

其複雜度

$O({\log }_{2}n)$

練習：

ZERO JUDGE b525 先別管這個了，你聽過turtlebee嗎？

---

質數判斷

若數

$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$

for (int i = 2; i <= sqrt(n); i++)
  if (n%i == 0) return false;

return true;

練習：

CODEFORCES 1033B Square Difference

Fermat's Primality Test

根據費馬小定理，

$$n\text{ is prime }\Rightarrow a^{n-1}\equiv 1\mathrm{mod}n\text{, where }n\nmid a$$

雖然根據邏輯規則以下逆命題不見得成立

$$n\text{ is prime }\Leftarrow a^{n-1}\equiv 1\mathrm{mod}n\text{, where }n\nmid a$$
但卻有很高的機率^[1]是成立的！

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 函數是快速冪演算法，則複雜度

$O(\log n)$

根據邏輯規則，如果回傳 false，那麼保證

$n$ 不是質數

$p\Rightarrow q$ 與

$\mathrm{\neg }p\Leftarrow \mathrm{\neg }q$ 是等價的

注意這是測試法，有時會遇到難判斷的數，例如

$n=561$ 對大部分的

$a$ 回傳 true

練習：

UVa OJ 10006 Carmichael Numbers

Miller Rabin Primality Test

Deus ex machina

繼續從 Fermat's Primality Test 發展：

除

$2$ 以外的質數都為奇數，故只關注奇數

$n$ 就好
則

$n-1$ 是個偶數，可將其寫成

$n-1=2^t\cdot d$

再注意到，若

$i<t$ 有

$a^{2^i\cdot d}\equiv \pm 1\mathrm{mod}n$ 則

$a^{2^t\cdot d}\equiv 1\mathrm{mod}n$

因為

$(\pm 1{)}^2$ 等於

$1$ 嘛

並且若

$n$ 是個質數，則

$$\begin{array}{rl}{x}^2\equiv 1\mathrm{mod}n& \Rightarrow n\mid x^2-1=(x+1)\cdot (x-1)\\ & \Rightarrow n\mid x+1\vee n\mid x-1\\ & \Rightarrow x\equiv -1\mathrm{mod}n\vee x\equiv 1\mathrm{mod}n\end{array}$$

如果

$n$ 是合數，第二列

$"\Rightarrow "$ 不一定會成立

有了上述性質，又能使質數測試成功機率大大提升

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)$ 完美判斷常用質數了！

練習：

ZERO JUDGE a007 判斷質數 ^[2]

---

質因數分解

根據唯一分解定理，任何合數都能被分解成一些質數的積

算術基本定理

$[2,\sqrt{n}]$ 試除法

與

$[2,\sqrt{n}]$ 判斷法原理相同

合數

$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}]$ 從小至大找數試除即可。

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
CODEFORCES 1114C Trailing Loves (or L'oeufs?)
LeetCode 952 Largest Component Size by Common Factor

---

質數篩檢

數質數可以有效安撫緊張的情緒

顧名思義，通常會想知道在一定的範圍內哪些數是質數

Sieve of Eratosthenes

其精神是將質數的倍數都設為非質數(合數)

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 溢位

練習：

UVa OJ 543 Goldbach’s Conjecture
UVa OJ 10140 Prime Distance

線性篩法

線性是指演算法時間

$T(n)=O(n)$

與 Sieve of Eratosthenes 相反：刪所有質數的倍數，就是刪所有數的質數倍
而刪質數倍的過程能簡單的判斷是否能及早收手，讓計算效率大幅提升！

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

存在性

數到

$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$

唯一性

任意

$n$ 只會被判定過恰好一次

令

$x_i=\frac{n}{p_i}$，

$i\in \{1,2,\cdots ,k\}$

$n$ 只會被數到

$x_1$ 時給判定，同樣用

$p_1$ 湊得

$n=p_1\cdot x_1$
對於

$x_i\ne x_1$， 在數到

$p_1$ 時

$x_i$ 就被

$p_1$ 給整除了，因而不會湊得

$p_i\cdot x_i$

$p_i$ 根據假設有

$p_i>p_1$

---

1. 此機率是根據測試的次數計算，每次測試會有不同的

   $a$ ↩︎

2. 這根本不是基礎題= = ↩︎