质数与GCD

发表于 分类于 阅读次数： Disqus： 本文字数： 2.5k 阅读时长 ≈ 9 分钟

简介

质数和最大公约数（GCD）是数论中最基础的部分之一。
本文将会讨论质数的判定与筛法以及欧几里得算法与扩展欧几里得算法。

质数

复习一下小学知识：
设\(a, b\)满足\(a \in \mathbb{N}^{\ast}, b \in \mathbb{N}\)。如果存在\(k \in \mathbb{N}\)使得\(b = k a\)，那么\(b\)是\(a\)的倍数（multiple），\(a\)是\(b\)的约数（divisor）或因数（factor）。这种关系记作\(a \mid b\)，读作\(a\)整除（divide）\(b\)。所有整数都整除\(0\)。
显然，任意大于\(1\)的正整数\(n\)都可以被\(1\)和\(n\)整除。如果除此之外\(n\)没有其它约数，那么\(n\)为质数（prime number）或素数，否则\(n\)为合数（composite number）。\(1\)既非质数，也非合数。

判定

暴力

要判定一个数是不是质数，只要枚举每个数看看是否可以整除就行了： 

1
2
3
4
5
6
7
8
bool 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
8
bool 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
11
bool 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
23
bool 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
12
int 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
16
int 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\)互质（coprime）。

欧几里得算法

如果我们已知两个数\(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\)的公约数相等，两者的最大公约数也必然相等。所以得到式子： \[\gcd(a, b) = \gcd(b, a \bmod b)\]

1
2
3
4
5
int gcd(int a, int b) {
    if (b == 0)
        return a;
    return gcd(b, a % 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
10
int 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
8
int 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
2
3
4
int inv(int a, int p, int &x, int &y) {
    exgcd(a, p, x, y);
    return (x % p + p) % p;
}

快速幂

这个方法需要运用费马小定理。费马小定理中：

若\(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
2
3
4
5
6
7
int inv(int a, int p) {
    int power = p - 2, res = 1;
    for (; power != 0; power >>= 1, a = (a * a) % p)
        if (power & 1)
            res = (a * res) % p;
    return res;
}

线性算法

对于求一连串数字的逆元，就只能用线性算法（\(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}\)。
然后我们就可以从前面的逆元推出后面的了： \[\operatorname{inv}_{i} = -\lfloor \dfrac{p}{i} \rfloor \operatorname{inv}_{p \bmod i} \bmod p\]

代码很短： 

1
2
3
inv[1] = 1;
for (int i = 2; i <= n; ++i)
    inv[i] = -(p / i) * inv[p % i] % p;
由于这种算法可能会算出负数，如果只求正整数，稍微修改一下代码即可： 
1
2
3
inv[1] = 1;
for (int i = 2; i <= n; ++i)
    inv[i] = (p - p / i) * inv[p % i] % p;

例题