简介
质数和GCD是数论中最基础的部分之一。
本篇将会讨论质数的判定与筛法，以及欧几里得算法与扩展欧几里得算法。
质数
复习一下小学知识：
设$a, b$满足$a \in \mathbb{N}^{\ast}, b \in \mathbb{N}$。若存在$k \in \mathbb{N}$使得$b = k a$，那么$b$是$a$的倍数，$a$是$b$的约数（也称因数）。这种关系记作$a \mid b$，读作$a$整除$b$。所有整数都整除$0$。
显然，任意大于$1$的正整数$n$都可以被$1$和$n$整除，若除此之外$n$没有其它约数，那么$n$为质数（也称素数），否则$n$为合数。$1$既非质数，也非合数。
判定
暴力
要判定一个数是不是质数，只要枚举每个数看看是否可以整除就行了：1
2
3
4
5
6
7
8bool isPrime(int n) {
if (n < 2)
return false;
for (int i = 2; i < n; ++i)
if (n % i == 0)
return false;
return true;
}
容易想到，如果$a$是$n$的约数，那么$\dfrac{n}{a}$也是$n$的约数。
因此，对于每一对$(a, \dfrac{n}{a})$，只要检验其中的一个就好了。每组数中的较小数都在区间$[1, \sqrt{n}]$中：1
2
3
4
5
6
7
8bool isPrime(int n) {
if (n < 2)
return false;
for (int i = 2; i <= sqrt(n); ++i)
if (n % i == 0)
return false;
return true;
}
Miller-Rabin
Miller-Rabin素性测试是一种更好的质数判定方法。
它对数$n$进行$k$轮测试的时间复杂度是$O(k \log^{3} n)$。
费马素性测试
费马小定理有：
若$p$为质数，对于任意整数$a$，$a^{p} \equiv a \pmod{p}$。
这个定理也可以表示成：
若$p$为质数，对于任意与$p$互质的整数$a$，$a^{p - 1} \equiv 1 \pmod{p}$。
我们可以在$[2, n - 1]$中随机选取整数$a$，检验是否$a^{n - 1} \equiv 1 \pmod{n}$：1
2
3
4
5
6
7
8
9
10
11bool fermat(int n) {
if (n <= 2)
return n == 2;
// k为测试次数，k越大结果越准确，但过大会影响效率
for (int i = 1; i <= k; ++i) {
int a = rand() % (n - 2) + 2;
if (quickPow(a, n - 1, n) != 1)
return false;
}
return true;
}
费马小定理的逆定理并不成立，换而言之，即使满足了$a^{n - 1} \equiv 1 \pmod{n}$，$n$也不一定是质数。但是上面的做法中随机选择$a$，大大降低了判断错误的概率。
不过仍有一类数会导致错误判断：
对于合数$n$，若对于所有与$a$互质的正整数$a$，$a^{n - 1} \equiv 1 \pmod{n}$都成立，则这个合数$n$称为卡迈克尔数（Carmichael），也称费马伪素数。
例如$561 = 3 \times 11 \times 17$就是一个卡迈克尔数。
二次探测定理
若$p$为奇质数，则$x^{2} \equiv 1 \pmod{p}$的解为$x = 1$或$x = p - 1$。
我们可以结合费马小定理和二次探测定理：
将$n - 1$分解为$n - 1 = u \times 2^{t}$，$u$为奇数。不断平方$u$，如果有非平凡平方根那么就不是质数：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23bool millerRabin(int n) {
if (n <= 2)
return n == 2;
int u = n - 1, t = 0;
while (u % 2 == 0) {
++t;
u /= 2;
}
// k为测试次数，k越大结果越准确，但过大会影响效率
for (int i = 1, j; i <= k; ++i) {
int a = rand() % (n - 2) + 2, b = pow(a, u, n);
if (b == 1 || b == n - 1)
continue;
for (j = 0; j < t; ++j) {
b = b * b % n;
if (b == n - 1)
break;
}
if (j >= t)
return false;
}
return true;
}
筛法
我们需要知道小于等于$n$有多少个质数。
最暴力的方法是对每个小于等于$n$的数进行判定，但是这样显然不能到达最优复杂度，于是考虑优化。
埃氏筛
如果$n$是合数，那么$n$的倍数也一定是合数。我们可以利用这个结论避免很多不必要的检测：1
2
3
4
5
6
7
8
9
10
11
12int eratosthenes(int n) {
int tot = 0, prime[MAXN];
bool vis[MAXN] = {};
for (int i = 2; i <= n; ++i) {
if (!vis[i]) {
prime[++tot] = i;
for (int j = 2 * i; j <= n; j += i)
vis[j] = true;
}
}
return tot;
}
以上为埃氏筛法，即埃拉托色尼筛法（Eratosthenes筛法），时间复杂度为$O(n \log \log n)$。
埃氏筛法仍然有优化空间，因为有的数被重复标记了。
欧拉筛
欧拉筛法（Euler筛法）是一种线性筛法。
欧拉筛中每个合数只被标记了一次，因此时间复杂度降到了$O(n)$：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16int euler(int n) {
int tot = 0, prime[MAX_SIZE];
bool vis[MAX_SIZE] = {};
for (int i = 2; i <= n; ++i) {
if (!vis[i])
prime[++tot] = i;
for (int j = 1; j <= tot; ++j) {
if (i * prime[j] > n)
break;
vis[i * prime[j]] = true;
if (i % prime[j] == 0)
break;
}
}
return tot;
}
例题
最大公约数
最大公约数（greatest common divisor）是指一组数中，每一个数都同时拥有的最大的约数，常缩写为GCD。
若$\gcd(a, b) = 1$，我们称$a$与$b$互质。
欧几里得算法
如果我们已知两个数$a, b$，如何求它们的最大公约数呢？
不妨设$a > b$，
若$b$为$a$的约数，那么$b$就是两者的最大公约数；
若$b$不能整除$a$，即$a = k b + c \ (c < b)$时：
显然$c = a \bmod b$，设$d \mid a, d \mid b$，则$c = a - k b, \dfrac{c}{d} = \dfrac{a}{d} - k \dfrac{b}{d}$。易知$\dfrac{c}{d}$为整数，即$d \mid c$。所以$d$是$b$和$a \bmod b$的公约数。
反过来，设$d \mid b, d \mid (a \bmod b)$，同样得出以下式子$\dfrac{a \bmod b}{d} = \dfrac{a}{d} - k \dfrac{b}{d}, \dfrac{a}{d} = \dfrac{a \bmod b}{d} + k \dfrac{b}{d}$。因为右边为整数，所以$\dfrac{a}{d}$也是整数，即$d \mid a$。因此$d$是$a$和$b$的公约数。
因$a, b$与$b, a \bmod b$的公约数相等，两者的最大公约数也必然相等。
所以得到式子：
1 | int gcd(int a, int b) { |
以上算法为欧几里得算法（Euclidean algorithm）。
扩展欧几里得算法
扩展欧几里得算法（extended Euclidean algorithm, EXGCD），常用于求$a x + b y = \gcd(a, b)$的一组可行解。
设$a x_{1} + b y_{1} = \gcd(a, b), b x_{2} + (a \bmod b) y_{2} = \gcd(b, a \bmod b)$，
由欧几里得定理可知：$\gcd(a, b) = \gcd(b, a \bmod b)$，
所以$a x_{1} + b y_{1} = b x_{2} + (a \bmod b) y_{2}$。
又因$a \bmod b = a - (\lfloor \dfrac{a}{b} \rfloor \times b)$，
所以$a x_{1} + b y_{1} = b x_{2} + (a - (\lfloor \dfrac{a}{b} \rfloor \times b)) y_{2}$。
$a x_{1} + b y_{1} = a y_{2} + b x_{2} - \lfloor \dfrac{a}{b} \rfloor \times b y_{2} = a y_{2} + b (x_{2} - \lfloor \dfrac{a}{b} \rfloor y_{2})$，
所以$x_{1} = y_{2}, y_{1} = x_{2} - \lfloor \dfrac{a}{b} \rfloor y_{2}$。
将$x_{2}, y_{2}$不断代入递归直到求出最大公约数为止，返回$x = 1, y = 0$求解，函数返回值为最大公约数：1
2
3
4
5
6
7
8
9
10int exgcd(int a, int b, int &x, int &y) {
if (b == 0) {
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;
}
线性同余方程
线性同余方程是最基础的同余方程，其未知数次数为一次，即形如$a x \equiv b \pmod{c}$的方程。
$a x \equiv b \pmod{c}$的本质其实就是$a x + c y = b$，根据裴蜀定理，当且仅当$\gcd(a, c) \mid b$时该方程存在整数解。
我们可以先用扩展欧几里得算法求出一组方程$a x_{0} + c y_{0} = \gcd(a, c)$的解$x_{0}, y_{0}$，然后两边同时除以$\gcd(a, c)$并乘上$b$，也就是$\dfrac{a b x_{0}}{\gcd(a, c)} + \dfrac{b c y_{0}}{\gcd(a, c)} = b$，我们就得到了方程的一组特解。
此外，若$x_{0}, y_{0}$是方程$a x + c y = b$的解，那么该方程的任意解可表示为$x = x_{0} + c t, y = y_{0} - a t$，对任意整数$t$都成立。
这样，我们就能得到所有解了。但是，题目往往都会让我们求最小正整数解，方法很简单：$t = \dfrac{c}{\gcd(a, c)}, x = (x \bmod t + t) \bmod t$。
代码：1
2
3
4
5
6
7
8int eq(int a, int b, int c, int &x, int &y) {
int d = exgcd(a, c, x, y);
if (b % d != 0)
return -1;
int k = b / d, t = c / d;
x *= k, y *= k;
return (x % t + t) % t;
}
乘法逆元
如果一个线性同余方程$a x \equiv 1 \pmod{p}$，那么$x$称为$a$模$p$意义下的逆元，记作$a^{-1}$。
扩展欧几里得
用扩展欧几里得的解法和解线性同余方程类似，就不多解释了。
1 | int inv(int a, int p, int &x, int &y) { |
快速幂
这个方法需要运用费马小定理。
费马小定理中：
若$p$为质数，对于任意与$p$互质的整数$a$，$a^{p - 1} \equiv 1 \pmod{p}$。
我们发现，右边刚好为$1$。
因为$a x \equiv 1 \pmod{p}$，
所以$a x \equiv a^{p - 1} \pmod{p}$，
所以$x \equiv a^{p - 2} \pmod{p}$。
然后我们就可以用快速幂求出$a^{p - 2} \bmod p$的值了，结果就是它的逆元。
1 | int inv(int a, int p) { |
线性算法
对于求一连串数字的逆元，就只能用线性算法（$O(n)$）。
首先，显然地，$1^{-1} \equiv 1 \pmod{p}$，
然后设$p = k i + r \ (1 < r < i < p)$，
将这个式子放到$\bmod p$意义下就会得到：$k i + r \equiv 0 \pmod{p}$。
两边同时乘$i^{-1}, r^{-1}$，得到$k r^{-1} + i^{-1} \equiv 0 \pmod{p}$，
$i^{-1} \equiv -k r^{-1} \pmod{p}$，
$i^{-1} \equiv -\lfloor \dfrac{p}{i} \rfloor (p \bmod i)^{-1} \pmod{p}$。
然后我们就可以从前面的逆元推出后面的了：
代码很短：1
2
3inv[1] = 1;
for (int i = 2; i <= n; ++i)
inv[i] = -(p / i) * inv[p % i] % p;
由于这种算法可能会算出负数，如果只求正整数，稍微修改一下代码即可：1
2
3inv[1] = 1;
for (int i = 2; i <= n; ++i)
inv[i] = (p - p / i) * inv[p % i] % p;
