入门级-数论-素（质）数判定、筛法、GCD、EXGCD等

# 素数判定

素数分布比较稀疏。

素数计数函数：小于或等于 $x$ 的素数的个数，用 $\pi(x)$ 表示。随着 $x$ 的增大，有这样的近似结果：$\pi(x) \sim \frac{x}{\ln(x)}$。

---

## 素数判定

如何用计算机来判断一个数是不是素数呢？

---

### 实现

暴力做法自然可以枚举从小到大的每个数看是否能整除

```cpp
// C++ Version
bool isPrime(a) {
  if (a < 2) return 0;
  for (int i = 2; i < a; ++i)
    if (a % i == 0) return 0;
  return 1;
}
```

---

这样做是十分稳妥了，但是真的有必要每个数都去判断吗？

很容易发现这样一个事实：如果 $x$ 是 $a$ 的约数，那么 $\frac{a}{x}$ 也是 $a$ 的约数。

这个结论告诉我们，对于每一对 $(x, \frac{a}{x} )$，只需要检验其中的一个就好了。为了方便起见，我们之考察每一对里面小的那个数。不难发现，所有这些较小数就是 $[1, \sqrt{a}]$ 这个区间里的数。

由于 $1$ 肯定是约数，所以不检验它。

---

```cpp
// C++ Version
bool isPrime(a) {
  if (a < 2) return 0;
  for (int i = 2; i * i <= a; ++i)
    if (a % i == 0) return 0;
  return 1;
}
```

---

# 素数筛法

---

## 1、埃氏筛法

根据质数的定义，能够被质数整除的数，肯定都是合数。

从小到大遍历每个数x，我们假设这个数是质数，把x的倍数都筛选掉，剩下的就是质数。

---

模拟筛法过程：

|i：|2|3|4|5|6|7|8|9|10|
|---|---|---|---|---|---|---|---|---|---|
|2:   |`2`|3|`4`|5|`6`|7|`8`|9|`10`|
|3:   |2|`3`|-|5|-|7|-|`9`|-|
|...|||||

---

此算法用到了双重循环，但由于内层循环的步长为i，时间复杂度可看成O(nloglogn)。

```
for(int i=2; i<=n; i++) {
    if (isPrime[i]) {//不是合数，也就是质数
      primes.push_back(i);
      for(int j=2*i; j<=n ; j=j+i)
       	 isPrime[j] = false;
    }
  }
```

---

## 2、欧拉筛法

在筛法模拟里会看到存在大量的重复筛，比如6即被2筛掉、也被3筛掉，
如果每个合数都只被一个数筛掉，那么这个算法显然就是O(n)的算法，即线性筛法。

`数学描述`：每个数都只被它的最小质因子筛掉。

---

`算法描述`：

* 对每个数a，将a之前的质数乘以a筛为合数。
* 如果某个质数是a的因子，则停止（后面的不要筛了）

对于某个数a，假设已经找到m个质数：p、p1、...pm <= a，   计划将 a 和 这些质数的乘积都筛掉。

---

但如果a存在一个质因子，假设为px，即 a = px * t
那么对于px到pm之间的任何质数p，都有 a * p = px * (t * p)
也就是说后面肯定会有一个比a大的数b = (t*p)，也会筛掉   b * px =  a * p

所以，这里就不需要筛这些质数（大于px的质数），可以直接跳过。

---

模拟筛法过程：

|i | 质数 | 筛
|- | - | -|
|2 | `2` | 4|
|3 | 2, `3` | 6, 9|
|4 | `2`,3 | 8|
|5| 2,3,`5` | 10, 15, 25|
|6| `2`,3,5 | 12|
|7|2,3,5,`7`|14,21,35,49|
|`8`|`2`,3,5,7|16|
|`9`|2,`3`,5,7|18,27|
|...|

---

注意上面的8和9，只需要筛到其质因子乘以该数为止。

---

由于每个数都只有其最小质因子乘以某个a筛掉，所以这个算法为O(n)。

```c++
for(int i=2; i<=n; i++) {
    if(isPrime[i])
      primes.push_back(i);
      for(int j=0; j<primes.size(); j++) {
        int p = primes[j];
        if(i * p > n) break;
        isPrime[i * p] = false;
        if(i%p == 0) break;
      }
  }
```

---

# GCD有关

## 定义

最大公约数即为 Greatest Common Divisor，常缩写为 gcd。

一组整数的公约数，是指同时是这组数中每一个数的约数的数。$\pm 1$ 是任意一组整数的公约数。

一组整数的最大公约数，是指所有公约数里面最大的一个。

---

那么如何求最大公约数呢？我们先考虑两个数的情况。

---

### 欧几里得算法

#### 过程

如果我们已知两个数 $a$ 和 $b$，如何求出二者的最大公约数呢？

不妨设 $a > b$

我们发现如果 $b$ 是 $a$ 的约数，那么 $b$ 就是二者的最大公约数。
下面讨论不能整除的情况，即 $a = b \times q + r$，其中 $r < b$。

我们通过证明可以得到 $\gcd(a,b)=\gcd(b,a \bmod b)$，过程如下：

---

###  "证明"
设 $a=bk+c$，显然有 $c=a \bmod b$。设 $d \mid a,~d \mid b$，则 $c=a-bk, \frac{c}{d}=\frac{a}{d}-\frac{b}{d}k$。
    
由右边的式子可知 $\frac{c}{d}$ 为整数，即 $d \mid c$，所以对于 $a,b$ 的公约数，它也会是 $a \bmod b$ 的公约数。

---

反过来也需要证明：
    
设 $d \mid b, \~d\mid (a \bmod b)$，我们还是可以像之前一样得到以下式子 $\frac{a\bmod b}{d}=\frac{a}{d}-\frac{b}{d}k,~\frac{a\bmod b}{d}+\frac{b}{d}k=\frac{a}{d}$。
    
因为左边式子显然为整数，所以 $\frac{a}{d}$ 也为整数，即 $d \mid a$，所以 $b,a\bmod b$ 的公约数也是 $a,b$ 的公约数。

---

既然两式公约数都是相同的，那么最大公约数也会相同。
    
所以得到式子 $\gcd(a,b)=\gcd(b,a\bmod b)$

既然得到了 $\gcd(a, b) = \gcd(b, r)$，这里两个数的大小是不会增大的，那么我们也就得到了关于两个数的最大公约数的一个递归求法。

---

#### 实现

```cpp
// C++ Version 1
int gcd(int a, int b) {
  if (b == 0) return a;
  return gcd(b, a % b);
}

// C++ Version 2
int gcd(int a, int b) { return b == 0 ? a : gcd(b, a % b); }
```

---

另外，对于 C++14，我们可以使用自带的 ` __gcd(a,b)` 函数来求最大公约数。而对于 C++ 17，我们可以使用 [`<numeric>`](https://en.cppreference.com/w/cpp/header/numeric) 头中的 [`std::gcd`](https://en.cppreference.com/w/cpp/numeric/gcd) 与 [`std::lcm`](https://en.cppreference.com/w/cpp/numeric/lcm) 来求最大公约数和最小公倍数。

---

递归至 `b==0`（即上一步的 `a%b==0`) 的情况再返回值即可。

上述算法被称作欧几里得算法（Euclidean algorithm）。

如果两个数 $a$ 和 $b$ 满足 $\gcd(a, b) = 1$，我们称 $a$ 和 $b$ 互质。

---

#### 性质

欧几里得算法的时间效率如何呢？下面我们证明，欧几里得算法的时间复杂度为 $O(\log n)$。
---

当我们求 $\gcd(a,b)$ 的时候，会遇到两种情况：
    
- $a < b$，这时候 $\gcd(a,b)=\gcd(b,a)$；
- $a \geq b$，这时候 $\gcd(a,b)=\gcd(b,a \bmod b)$，而对 $a$ 取模会让 $a$ 至少折半。这意味着这一过程最多发生 $O(\log n)$ 次。
    
第一种情况发生后一定会发生第二种情况，因此第一种情况的发生次数一定 **不多于** 第二种情况的发生次数。
    
从而我们最多递归 $O(\log n)$ 次就可以得出结果。

事实上，假如我们试着用欧几里得算法去求 [斐波那契数列](../combinatorics/fibonacci.md) 相邻两项的最大公约数，会让该算法达到最坏复杂度。

---

### 多个数的最大公约数

那怎么求多个数的最大公约数呢？显然答案一定是每个数的约数，那么也一定是每相邻两个数的约数。我们采用归纳法，可以证明，每次取出两个数求出答案后再放回去，不会对所需要的答案造成影响。

---

## 最小公倍数

接下来我们介绍如何求解最小公倍数（Least Common Multiple, LCM）。

---

### 定义

一组整数的公倍数，是指同时是这组数中每一个数的倍数的数。0 是任意一组整数的公倍数。

一组整数的最小公倍数，是指所有正的公倍数里面，最小的一个数。

---

### 两个数

设 $a = p_1^{k_{a_1}}p_2^{k_{a_2}} \cdots p_s^{k_{a_s}}$，$b = p_1^{k_{b_1}}p_2^{k_{b_2}} \cdots p_s^{k_{b_s}}$

我们发现，对于 $a$ 和 $b$ 的情况，二者的最大公约数等于

$p_1^{\min(k_{a_1}, k_{b_1})}p_2^{\min(k_{a_2}, k_{b_2})} \cdots p_s^{\min(k_{a_s}, k_{b_s})}$

---

最小公倍数等于

$p_1^{\max(k_{a_1}, k_{b_1})}p_2^{\max(k_{a_2}, k_{b_2})} \cdots p_s^{\max(k_{a_s}, k_{b_s})}$

由于 $k_a + k_b = \max(k_a, k_b) + \min(k_a, k_b)$

所以得到结论是 $\gcd(a, b) \times \operatorname{lcm}(a, b) = a \times b$

要求两个数的最小公倍数，先求出最大公约数即可。

---

### 多个数

可以发现，当我们求出两个数的 $\gcd$ 时，求最小公倍数是 $O(1)$ 的复杂度。那么对于多个数，我们其实没有必要求一个共同的最大公约数再去处理，最直接的方法就是，当我们算出两个数的 $\gcd$，或许在求多个数的 $\gcd$ 时候，我们将它放入序列对后面的数继续求解，那么，我们转换一下，直接将最小公倍数放入序列即可。

---

## 扩展欧几里得算法

扩展欧几里得算法（Extended Euclidean algorithm, EXGCD），常用于求 $ax+by=\gcd(a,b)$ 的一组可行解。

---

### 过程

设

$ax_1+by_1=\gcd(a,b)$

$bx_2+(a\bmod b)y_2=\gcd(b,a\bmod b)$

由欧几里得定理可知：$\gcd(a,b)=\gcd(b,a\bmod b)$

所以 $ax_1+by_1=bx_2+(a\bmod b)y_2$

又因为 $a\bmod b=a-(\lfloor\frac{a}{b}\rfloor\times b)$

----

所以 $ax_1+by_1=bx_2+(a-(\lfloor\frac{a}{b}\rfloor\times b))y_2$

$ax_1+by_1=ay_2+bx_2-\lfloor\frac{a}{b}\rfloor\times by_2=ay_2+b(x_2-\lfloor\frac{a}{b}\rfloor y_2)$

因为 $a=a,b=b$，所以 $x_1=y_2,y_1=x_2-\lfloor\frac{a}{b}\rfloor y_2$
将 $x_2,y_2$ 不断代入递归求解直至 $\gcd$（最大公约数，下同）为 `0` 递归 `x=1,y=0` 回去求解。

---

### 实现

```cpp
// C++ Version
int Exgcd(int a, int b, int &x, int &y) {
  if (!b) {
    x = 1;
    y = 0;
    return a;
  }
  int d = Exgcd(b, a % b, x, y);
  int t = x;
  x = y;
  y = t - (a / b) * y;
  return d;
}
```

函数返回的值为 $\gcd$，在这个过程中计算 $x,y$ 即可。

---

### 值域分析

$ax+by=\gcd(a,b)$ 的解有无数个，显然其中有的解会爆 long long。  
万幸的是，若 $b\not= 0$，扩展欧几里得算法求出的可行解必有 $|x|\le b,|y|\le a$。  
下面给出这一性质的证明。

"证明"
-   $\gcd(a,b)=b$ 时，$a\bmod b=0$，必在下一层终止递归。  
        得到 $x_1=0,y_1=1$，显然 $a,b\ge 1\ge |x_1|,|y_1|$。
-   $\gcd(a,b)\not= b$ 时，设 $|x_2|\le (a\bmod b),|y_2|\le b$。  
因为 $x_1=y_2,y_1=x_2-{\left\lfloor\dfrac{a}{b}\right\rfloor}y_2$   
所以 $|x_1|=|y_2|\le b,|y_1|\le|x_2|+|{\left\lfloor\dfrac{a}{b}\right\rfloor}y_2|\le (a\bmod b)+{\left\lfloor\dfrac{a}{b}\right\rfloor}|y_2|$  
$\le a-{\left\lfloor\dfrac{a}{b}\right\rfloor}b+{\left\lfloor\dfrac{a}{b}\right\rfloor}|y_2|\le a-{\left\lfloor\dfrac{a}{b}\right\rfloor}(b-|y_2|)$   
$a\bmod b=a-{\left\lfloor\dfrac{a}{b}\right\rfloor}b\le a-{\left\lfloor\dfrac{a}{b}\right\rfloor}(b-|y_2|)\le a$   
因此 $|x_1|\le b,|y_1|\le a$ 成立。

---

### 迭代法编写扩展欧几里得算法

因为迭代的方法避免了递归，所以代码运行速度将比递归代码快一点。

```cpp
int gcd(int a, int b, int& x, int& y) {
  x = 1, y = 0;
  int x1 = 0, y1 = 1, a1 = a, b1 = b;
  while (b1) {
    int q = a1 / b1;
    tie(x, x1) = make_tuple(x1, x - q * x1);
    tie(y, y1) = make_tuple(y1, y - q * y1);
    tie(a1, b1) = make_tuple(b1, a1 - q * b1);
  }
  return a1;
}
```

如果你仔细观察 $a_1$ 和 $b_1$，你会发现，他们在迭代版本的欧几里德算法中取值完全相同，并且以下公式无论何时（在 while 循环之前和每次迭代结束时）都是成立的：$x \cdot a +y \cdot b =a_1$ 和 $x_1 \cdot a +y_1 \cdot b= b_1$。因此，该算法肯定能正确计算出 $\gcd$。

最后我们知道 $a_1$ 就是要求的 $\gcd$，有 $x \cdot a +y \cdot b =g$。

---

#### 矩阵的解释

对于正整数 $a$ 和 $b$ 的一次辗转相除即 $\gcd(a,b)=\gcd(b,a\bmod b)$ 使用矩阵表示如

$$
\begin{bmatrix}
b\\\\
a\bmod b
\end{bmatrix}
\=
\begin{bmatrix}
0\&1\\\\
1\&-\lfloor a/b\rfloor
\end{bmatrix}
\begin{bmatrix}
\&a\\\\
\&b
\end{bmatrix}
$$

其中向下取整符号 $\lfloor c\rfloor$ 表示不大于 $c$ 的最大整数。我们定义变换
$\begin{bmatrix}a\\\\b\end{bmatrix}\mapsto \begin{bmatrix}0&1\\\\1&-\lfloor a/b\rfloor\end{bmatrix}\begin{bmatrix}a\\\\b\end{bmatrix}$。

---

易发现欧几里得算法即不停应用该变换，有

$$
\begin{bmatrix}
\gcd(a,b)\\\\0
\end{bmatrix}
\=
\left(
\cdots 
\begin{bmatrix}
0&1\\\\1&-\lfloor a/b\rfloor
\end{bmatrix}
\begin{bmatrix}
1&0\\\\0&1
\end{bmatrix}
\right)
\begin{bmatrix}
a\\\\b
\end{bmatrix}
$$

---

令

$$
\begin{bmatrix}
x_1&x_2\\\\x_3&x_4
\end{bmatrix}
\=
\cdots 
\begin{bmatrix}
0&1\\\\1&-\lfloor a/b\rfloor
\end{bmatrix}
\begin{bmatrix}
1&0\\\\0&1
\end{bmatrix}
$$

那么

$$
\begin{bmatrix}
\gcd(a,b)\\\\0
\end{bmatrix}
\=
\begin{bmatrix}
x_1&x_2\\\\x_3&x_4
\end{bmatrix}
\begin{bmatrix}
a\\\\b
\end{bmatrix}
$$

---

满足 $a\cdot x_1+b\cdot x_2=\gcd(a,b)$ 即扩展欧几里得算法，注意在最后乘了一个单位矩阵不会影响结果，提示我们可以在开始时维护一个 $2\times 2$ 的单位矩阵编写更简洁的迭代方法如

```cpp
int exgcd(int a, int b, int &x, int &y) {
  int x1 = 1, x2 = 0, x3 = 0, x4 = 1;
  while (b != 0) {
    int c = a / b;
    std::tie(x1, x2, x3, x4, a, b) =
        std::make_tuple(x3, x4, x1 - x3 * c, x2 - x4 * c, b, a - b * c);
  }
  x = x1, y = x2;
  return a;
}
```

这种表述相较于递归更简单。