ACM
February 3, 2021

2021.2.3

数论

素数

简介

一个大于 1 的自然数，除了 1 和它自身外，不能被其他自然数整除的数叫做素数；否则称为合数

判断方法

暴力 · $O(\frac{n}{2})$

1
2
3
4
5
6
7
bool isPrime(a)
{
	if (a < 2) return 0;
	for (int i = 2; i * i < a; ++i)
		if (!(a % i)) return 0;
	return 1;
}

Miller-Rabin 素性测试 · $ O(k \ log^{3}n)$

Fermat 素性测试

我们可以根据费马小定理得出一种检验素数的思路：

它的基本思想是不断地选取在 [2, n - 1] 中的基 a，并检验是否每次都有 $a^{n - 1} \ \equiv \ 1 \ (mod \ n)$

1
2
3
4
5
6
7
8
9
10
11
12
bool millerRabin(int n)
{
    if (n < 3) return n == 2;
    // test_time 为测试次数,建议设为不小于 8
    // 的整数以保证正确率,但也不宜过大,否则会影响效率
    for (int i = 1; i <= test_time; ++i)
    {
        int a = rand() % (n - 2) + 2;
        if (quickPow(a, n - 1, n) != 1) return 0;
    }
    return 1;
}

很遗憾，费马小定理的逆定理并不成立，换言之，满足了 $a^{n - 1} \ \equiv \ 1 \ (mod \ n)$，a 也不一定是素数。

卡迈克尔数

上面的做法中随机地选择 a，很大程度地降低了犯错的概率。但是仍有一类数，上面的做法并不能准确地判断。

对于合数 n，如果对于所有正整数 a，a 和 n 互素，都有同余式 $a^{n - 1} \ \equiv \ 1 \ (mod \ n)$ 成立，则合数 n 为卡迈克尔数（Carmichael Number），又称为费马伪素数。

比如，561 = 3 x 11 x 17 就是一个卡迈克尔数。

而且我们知道，若 n 为卡迈克尔数，则 m = 2n - 1 也是一个卡迈克尔数，从而卡迈克尔数的个数是无穷的。

二次探测定理

如果 p 是奇素数，则 $x^{2} \ \equiv \ 1 \ (mod \ p)$ 的解为 $x \ \equiv \ 1 \ (mod \ p)$ 或者 $x \ \equiv \ p-1 \ (mod \ p)$ 。

要证明该定理，只需将上面的方程移项，再使用平方差公式，得到 $(x + 1)(x - 1) \ \equiv \ 0 \ (mod \ p)$，即可得出上面的结论。

实现

根据卡迈克尔数的性质，可知其一定不是 pe。

不妨将费马小定理和二次探测定理结合起来使用：

将 $a^{n - 1} \ \equiv \ 1 \ (mod \ n)$ 中的指数 n - 1 分解为 $n \ - \ 1 \ = \ u \ \times \ 2^{t}$，在每轮测试中对随机出来的 a 先求出 $a^{u} \ (mod \ n)$，之后对这个值执行最多 t 次平方操作，若发现非平凡平方根时即可判断出其不是素数，否则通过此轮测试。

模板
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
bool millerRabbin(int n)
{
    if (n < 3) return n == 2;
    int a = n - 1, b = 0;
    while (!(a % 2))
    {
        a /= 2;
        ++b;
    }
    // test_time 为测试次数,建议设为不小于 8 的整数以保证正确率,但也不宜过大,否则会影响效率
    for (int i = 1, j; i <= test_time; ++i)
    {
        int x = rand() % (n - 2) + 2, v = quickPow(x, a, n);
        if (v == 1 || v == n - 1) continue;
        for (j = 0; j < b; ++j)
        {
            v = (long long)v * v % n;
            if (v == n - 1) break;
        }
        if (j >= b) return 0;
    }
    return 1;
}

最大公约数 · GCD

简介

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

在素数一节中，我们已经介绍了约数的概念。

一组数的公约数，是指同时是这组数中每一个数的约数的数。而最大公约数，则是指所有公约数里面最大的一个。

求法

欧几里得算法 · $O(\log n)$

1
2
3
4
5
int gcd(int a, int b)
{
	if (!b) return a;
	return gcd(b, a % b);
}

多个数的最大公约数

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
int a[MAXN];
int multi_gcd(int n)
{
    int tmp1, tmp2;
    tmp1 = a[1];
    for (int i = 2; i <= n; i++)
    {
        tmp2 = a[i];
        a[i] = gcd(tmp1, tmp2);
        tmp1 = a[i];
    }

    return a[n];
}

最小公倍数 · LCM

1
2
3
4
int lcm(int a, int b)
{
	return a / gcd(a, b) * b;	// 先除再乘避免溢出
}

多个数的最小公倍数

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
int a[MAXN];
int multi_lcm(int n)
{
    int tmp1, tmp2;
    tmp1 = a[1];
    for (int i = 2; i <= n; i++)
    {
        tmp2 = a[i];
        a[i] = lcm(tmp1, tmp2);
        tmp1 = a[i];
    }

    return a[n];
}

欧拉函数

定义

欧拉函数（Euler’s totient function），即 $\phi(n)$，表示的是小于等于 n 和 n 互质的数的个数。

比如说 $\phi(1) = 1$。

当 n 是质数的时候，显然有 $\phi(n) = n - 1$。

性质

求欧拉函数值

1
2
3
4
5
6
7
8
9
10
11
12
13
int euler_phi(int n)
{
    int m = int(sqrt(n + 0.5));
    int ans = n;
    for (int i = 2; i <= m; i++)
        if (n % i == 0)
        {
            ans = ans / i * (i - 1);
            while (n % i == 0) n /= i;
        }
    if (n > 1) ans = ans / n * (n - 1);
    return ans;
}

如果将上面的程序改成如下形式，会提升一点效率：

1
2
3
4
5
6
7
8
9
10
11
12
int euler_phi(int n)
{
    int ans = n;
    for (int i = 2; i * i <= n; i++)
        if (n % i == 0)
        {
            ans = ans / i * (i - 1);
            while (n % i == 0) n /= i;
        }
    if (n > 1) ans = ans / n * (n - 1);
    return ans;
}

欧拉定理

若 $gcd(a, \ m) \ = \ 1$ 则 $a^{\phi(m)} \equiv 1 (mod m)$

扩展欧拉定理

$$
a^{b} \equiv \left{
\begin{aligned}
a^{b \ mod \ \phi(p)} & , & gcd(a, \ p) \ = \ 1 \
a^{b} & , & gcd(a, \ p) \ \neq \ 1 & , \ b \ < \ \phi(p) \ (mod \ p)\
a^{b \ mod \ \phi(p)+\phi(p)} & , & gcd(a, \ b) \ \neq \ 1 & , \ b \ \geq \ \phi(p)
\end{aligned}
\right.
$$

筛法

素数筛法

如果我们想要知道小于等于 n 有多少个素数呢？

一个自然的想法是对于小于等于 n 的每个数进行一次质数检验。这种暴力的做法显然不能达到最优复杂度。

埃拉托斯特尼（Eratosthenes）筛法 · $O(n\log\log n)$

如果 x 是合数，那么 x 的倍数也一定是合数。利用这个结论，我们可以避免很多次不必要的检测。

如果我们从小到大考虑每个数，然后同时把当前这个数的所有（比自己大的）倍数记为合数，那么运行结束的时候没有被标记的数就是素数了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
int Eratosthenes(int n)
{
    int p = 0;
    for (int i = 0; i <= n; ++i) is_prime[i] = 1;
    is_prime[0] = is_prime[1] = 0;
    for (int i = 2; i <= n; ++i)
    {
        if (is_prime[i])
        {
            prime[p++] = i;  // prime[p]是i,后置自增运算代表当前素数数量
            if ((long long)i * i <= n)
                for (int j = i * i; j <= n; j += i)
                    // 因为从 2 到 i - 1 的倍数我们之前筛过了，这里直接从 i
                    // 的倍数开始，提高了运行速度
                    is_prime[j] = 0;  // 是i的倍数的均不是素数
        }
    }
    return p;
}
筛至平方根 && 只筛奇数

显然，要找到直到 n 为止的所有素数，仅对不超过 $\sqrt{n}$ 的素数进行筛选就足够了。

因为除 2 以外的偶数都是合数，所以我们可以直接跳过它们，只用关心奇数就好。

1
2
3
4
5
6
7
8
9
10
11
int n;
vector<char> is_prime(n + 1, true);
is_prime[0] = is_prime[1] = false;
for (int i = 2; i * i <= n; i++)
{
    if (!(i % 2)) continue;
    if (is_prime[i])
    {
        for (int j = i * i; j <= n; j += i) is_prime[j] = false;
    }
}

线性筛法

埃氏筛法仍有优化空间，它会将一个合数重复多次标记。有没有什么办法省掉无意义的步骤呢？答案是肯定的。

如果能让每个合数都只被标记一次，那么时间复杂度就可以降到 $O(n)$ 了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
void init()
{
    phi[1] = 1;
    for (int i = 2; i < MAXN; ++i)
    {
        if (!vis[i])
        {
            phi[i] = i - 1;
            pri[cnt++] = i;
        }
        for (int j = 0; j < cnt; ++j)
        {
            if (1ll * i * pri[j] >= MAXN) break;
            vis[i * pri[j]] = 1;
            if (i % pri[j])
            {
                phi[i * pri[j]] = phi[i] * (pri[j] - 1);
            }
            else
            {
                // i % pri[j] == 0
                // 换言之，i 之前被 pri[j] 筛过了
                // 由于 pri 里面质数是从小到大的，所以 i 乘上其他的质数的结果一定也是
                // pri[j] 的倍数 它们都被筛过了，就不需要再筛了，所以这里直接 break
                // 掉就好了
                phi[i * pri[j]] = phi[i] * pri[j];
                break;
            }
        }
    }
}

筛法求欧拉函数

1
2
3
4
5
6
7
8
9
10
11
12
void phi_table(int n, int* phi)
{
    for (int i = 2; i <= n; i++) phi[i] = 0;
    phi[1] = 1;
    for (int i = 2; i <= n; i++)
        if (!phi[i])
            for (int j = i; j <= n; j += i)
            {
                if (!phi[j]) phi[j] = j;
                phi[j] = phi[j] / i * (i - 1);
            }
}

筛法求约数个数

d[i]: i 的约数个数

num[i]: i 的最小质因子出现次数

定理：若 $n \ = \ \prod_{i=1}^{m} p_{i}^{c_{i}}$ 则 $d_{i} \ = \ \prod_{i=1}^{m}c_{i} \ + \ 1$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
void pre()
{
    d[1] = 1;
    for (int i = 2; i <= n; ++i)
    {
        if (!v[i]) v[i] = 1, p[++tot] = i, d[i] = 2, num[i] = 1;
        for (int j = 1; j <= tot && i <= n / p[j]; ++j)
        {
            v[p[j] * i] = 1;
            if (i % p[j] == 0)
            {
                num[i * p[j]] = num[i] + 1;
                d[i * p[j]] = d[i] / num[i * p[j]] * (num[i * p[j]] + 1);
                break;
            }
            else
            {
                num[i * p[j]] = 1;
                d[i * p[j]] = d[i] * 2;
            }
        }
    }
}

筛法求约数和

f[i]: i 的约数和

g[i]: i 的最小质因子的 $p \ + \ p^{1} \ + \ p^{2} \ + \ \cdots \ + \ p^{k}$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
void pre()
{
    g[1] = f[1] = 1;
    for (int i = 2; i <= n; ++i)
    {
        if (!v[i]) v[i] = 1, p[++tot] = i, g[i] = i + 1, f[i] = i + 1;
        for (int j = 1; j <= tot && i <= n / p[j]; ++j)
        {
            v[p[j] * i] = 1;
            if (i % p[j] == 0)
            {
                g[i * p[j]] = g[i] * p[j] + 1;
                f[i * p[j]] = f[i] / g[i] * g[i * p[j]];
                break;
            }
            else
            {
                f[i * p[j]] = f[i] * f[p[j]];
                g[i * p[j]] = 1 + p[j];
            }
        }
    }
    for (int i = 1; i <= n; ++i) f[i] = (f[i - 1] + f[i]) % Mod;
}

乘法逆元

简介

如果一个线性同余方程 $ax \ \equiv \ 1 \ (mod \ b)$ ，则 x 称为 a mod b 的逆元，记作 $a^{-1}$。

如何求逆元

扩展欧几里得法

1
2
3
4
5
6
7
8
9
10
11
void exgcd(int a, int b, int& x, int& y)
{
    if (!b)
    {
        x = 1;
        y = 0;
        return;
    }
    exgcd(b, a % b, y, x);
    y -= a / b * x;
}

快速幂法

因为 $ax \ \equiv \ 1 \ (mod \ b)$；

所以 $ax \ \equiv \ a^{b-1} \ (mod \ b)$（根据费马小定理）；

所以 $x \ \equiv \ a^{b-2} \ (mod \ b)$。

然后我们就可以用快速幂来求了。

1
2
3
4
5
6
7
8
9
10
11
inline int qpow(long long a, int b)
{
    int ans = 1;
    a = (a % p + p) % p;
    for (; b; b >>= 1)
    {
        if (b & 1) ans = (a * ans) % p;
        a = (a * a) % p;
    }
    return ans;
}

线性求逆元

1
2
3
4
5
inv[1] = 1;
for (int i = 2; i <= n; ++i)
{
	inv[i] = (long long)(p - p / i) * inv[p % i] % p;
}

线性求任意 n 个数的逆元

1
2
3
4
5
6
s[0] = 1;
for (int i = 1; i <= n; ++i) s[i] = s[i - 1] * a[i] % p;
sv[n] = qpow(s[n], p - 2);
// 当然这里也可以用 exgcd 来求逆元,视个人喜好而定.
for (int i = n; i >= 1; --i) sv[i - 1] = sv[i] * a[i] % p;
for (int i = 1; i <= n; ++i) inv[i] = sv[i] * s[i - 1] % p;

中国剩余定理

简介

中国剩余定理 (Chinese Remainder Theorem, CRT) 可求解如下形式的一元线性同余方程组（其中 $n_{1}, \ n_{2}, \ \cdots, \ n_{k}$ 两两互质）：
$$
\left {
\begin{array}{c}
x \ &\equiv \ &a_{1} \ (mod \ n_{1}) \
x \ &\equiv \ &a_{2} \ (mod \ n_{2}) \
&\cdots\ \
x \ &\equiv \ &a_{k} \ (mod \ n_{k})
\end{array}
\right .
$$

算法流程

  1. 计算所有模数的积 n；
  2. 对于第 i 个方程：
    1. 计算 $m_{i} \ = \ \frac{n}{n_{i}}$；
    2. 计算 $m_{i}$ 在模 $n_{i}$ 意义下的逆元 $m_{i}^{-1}$
    3. 计算 $c_{i} \ = \ m_{i}m_{i}^{-1}$（不要对 $n_{i}$ 取模）。
  3. 方程组的唯一解为：$a \ = \ \sum_{i=1}^{k}a_{i}c_{i} \ (mod \ n)$。

About this Post

This post is written by OwlllOvO, licensed under CC BY-NC 4.0.

#C++#数论