数论初步

数论主要研究整数的性质，被数学王子高斯称为“数学上的皇冠”。初等数论是用初等方法研究的数论，它的研究方法本质上说，就是利用整数环的整除性质，主要包括整除理论、同余理论、连分数理论。在 OI 中，几乎跟数学有关的题目都绕不开初等数论，本文将引导你学习初等数论的基础知识。

OI 中的数学并没有数竞中的那么难，因为正确证明不再是首要任务（尽管证明可以帮助你理解），应用才是王道。

在阅读本文前，请保证您学过小学数学（雾

如果你想以数学竞赛的模式学习数论，请参阅潘承洞的《初等数论》，并提前学习高等数学基本内容：

所以还是以 OI 的模式学习数论吧(￣ ‘i ￣;)。

注意，OI 的模式学习数论，是指证明过程不会那么详细，会显得“感性理解”，甚至有些直接跳过证明，但足以让你在考场上回忆起这些性质并用它们解题（考场上想的做法往往可以用对拍或其他手段来验证，感性证明“好像是对的”就足够了，重要的还是大量做题）。严谨甚至多种的证明过程请参考《初等数论》《具体数学》等书籍——这未必不是好事，虽然会花费时间，但收获也是巨大的。

希望在阅读时您能拿起纸笔和键盘，适当的演算和完成代码实现。

基础知识

接下来若没有特殊说明，我们的所有数都满足 $x \in \mathbb{Z}$。

模运算

大家都应该知道取模 $\text{mod}$。

随时取模性质

就是说 $\sum (x_i \bmod p) = (\sum x_i)\bmod p$，$\prod (x_i \bmod p) = (\prod x_i)\bmod p$。用人话说，运算的最终结果要取模，那么运算过程中的加法、减法和乘法都可以随时取模（但是除法不行！需要用下文介绍的逆元来解决）。

这东西可以证明，下文会出现（

C++ 的取模

在 C++ 中，用 % 来表示取模。

这里问个问题，$-5\bmod 3$ 是多少？

有的人会说是 $-2$，实际上是错误的。

我们不知道 $-5\bmod 3$ 是多少，但我们知道它等于 $1\bmod 3$，所以它的答案是 $1$。

从 C99 和 C++11 起，规定商向零取整（舍弃小数部分），取模的符号即与被除数相同。从此以下运算结果保证为真：

5 % 3 == 2;
5 % -3 == 2;
-5 % 3 == -2;
-5 % -3 == -2;

这也就意味着，如果存在负数，那么必须这样实现取模：(a % b + b) % b。

整除与约数

若 $a=bk$，则 $b$ 整除 $a$，记作 $b\mid a$。不整除记作 $b\nmid a$。
也称 $b$ 是 $a$ 的约数（因数），$a$ 是 $b$ 的倍数。

整除有一些性质：

• $a\mid b \iff -a\mid b \iff a\mid -b \iff |a| \mid |b|$；
• 如果 $a\mid b, a\mid c$，则可得 $a\mid b+c$；
• 如果 $a\mid b, b\mid c$，则 $a\mid c$，也就是整除的传递性；
• 若 $m\ne 0$，则 $a\mid b \iff ma\mid mb$；

证明如下：

还有：

• $a\mid b$ 且 $b\mid a \Longrightarrow b=\pm a$，这是因为这两个式子代表 $b=ap,a=bq\Rightarrow ab=abpq\Rightarrow pq=1\Rightarrow p=\pm 1$；
• 设 $b\neq 0$，则 $a\mid b \Longrightarrow |a|\le |b|$；
• 若 $a\neq 0,b=qa+c$，那么 $a\mid b \iff a\mid c$。

核心的证明思路就是将 $a\mid b$ 写作 $b=pa$。

求 N 的正约数集合

若 $a$ 是 $n$ 的约数，则 $\cfrac{n}{a}$ 也是 $n$ 的约数，且 $a$ 和 $\cfrac{n}{a}$ 之中必有一个 $\leqslant \sqrt{n}$。
枚举 $1\sim \sqrt{n}$ 的所有数，判断是不是 $n$ 的约数即可，代码如下：

int factor[1605];
inline int divisor(int n) // 返回约数个数，存在 factor 中
{
    int m = 0, t = min(n, sqrt(n) + 1);
    for (int i = 1; i <= t; ++i)
        if (n % i == 0)
        {
            factor[++m] = i;
            if (i != n / i) factor[++m] = n / i;
        }
    //sort(factor + 1, factor + m + 1);
    //此步有可能是必要的，比如在某些 DP 中，转移要有顺序。
    return m;
}

以上方法称为试除法，通过试除法可以证明，一个整数 $n$ 的约数个数（记作 $d(n)$）上界为 $2\sqrt{n}$（实际上远远不到，经过试验能得到下表）。

范围	$d(n)$ 最大的 $n$	$\max\{d(n)\}$
$\leqslant 10^4$	数据丢失	$64$
$\leqslant 10^5$	$83160$	$128$
$\leqslant 10^6$	$720720$	$240$
$\leqslant 10^7$	$8648640$	$448$
$\leqslant 10^8$	$73513440$	$768$
$\le 10^9$		$1344$
$\le 10^{12}$		$6720$
$\le 10^{15}$		$26880$
64 位整数	有这种数据吗	$103680$

“Pollard Rho” 是一个比“试除法”更高效的算法。难度大大超出了本文的范围，感兴趣的读者可以自行了解。

1~n 每个数的正约数集合

求 $1\sim n$ 每个数的正约数集合可以采用倍数法。对于每个数 $d$，$d,2d,3d,\cdots,\lfloor \cfrac{n}{d}\rfloor\times d$ 的约数就有 $d$。代码如下：

vector <int> factor[1000005];
for (int i = 1; i <= n; ++i)
    for (int j = i; j <= n; j += i)
        factor[j].push_back(i);

时间复杂度为 $O(\sum\limits_{i=1}^{n}\cfrac{n}{i})=O(n\log n)$，这一等式称之为调和级数。

GCD 与 LCM

若 $c$ 同时是 $a$ 和 $b$ 的约数，则称 $c$ 是 $a$ 和 $b$ 的公约数。
$a$ 和 $b$ 的最大公约数（Greatest Common Divisor，GCD）记为 $\gcd(a,b)$ ，或者在没有歧义的条件下记作 $(a,b)$。

若 $c$ 同时是 $a$ 和 $b$ 的倍数，则称 $c$ 是 $a$ 和 $b$ 的公倍数。
$a$ 和 $b$ 的最小公倍数（Least Common Multiple，LCM）记为 $\text{lcm}(a,b)$ ，或者在没有歧义的条件下记作 $[a,b]$。

GCD 和 LCM 有一些显然的性质，如下：

• $(a,a)=(0,a)=a,$
• $\forall a\mid b, (a,b)=a,$
• $(a,b)=(a,a+b)=(a,ka+b) \Rightarrow (a,b)=(a,b\bmod a),$
• $(ka,kb)=k\times(a,b),$
• $(a,b,c)=((a,b),c),$
• $[a,b,c]=[[a,b],c],$
• $[a,b]=a\times b \div (a,b),$
• 如果 $d$ 是 $a,b$ 的公约数，那么 $d\mid (a,b).$

值得注意的是，若 $(a,b)=1$，则称 $a$ 和 $b$ 互质（互素），记作 $a\perp b$^[1]。

还有一些比较难证的家伙，请关注一下：

• 如果 $(n,x)=1$，那么 $(n,n-x)=1$，
• 若 $a\perp m$，则 $(m,ab)=(m,b)$，
• 若 $a\perp m,m\mid ab$，则 $m\mid b$，
• $[a,b](a,b)=ab$。

证明：

完整的最大公约数理论及其不同的证明体系请参考《初等数论》。

求 GCD 与 LCM

求 GCD 有两种方法：辗转相除和更相减损。

辗转相除法

刚才我们知道，$(a,b)=(a,b\bmod a)$，但注意代码实现不能这么写！
比如给个 $(2,1)$，它就会一直是 $(2,1)$。

所以要这么写：

int gcd(int x, int y)
{
    if (y == 0) return x;
    return gcd(y, x % y);
}

看见了吗？反过来就行了。

那么它的复杂度是多少？

• 若 $a<b$，这时候 $(a,b)=(b,a)$；
• 若 $a\geq b$，这时候 $(a,b)=(b,a \bmod b)$，而对 $a$ 取模会让 $a$ 至少折半。这意味着这一过程最多进行 $\mathcal{O}(\log n)$ 次。

而不可能连续两次发生第一种情况。所以复杂度是 $\mathcal{O}(\log(a+b))$，但实际上远远达不到这个上界。当 $a=F_x,b=F_{x-1}$ 时（$F$ 指斐波那契数列），会跑满 $\mathcal{O}(\log(a+b))$，因为模运算会被卡成减法，也就成了更相减损法。

更相减损法

还是根据性质，有：

int gcd(int x, int y)
{
    if (y == 0) return x;
    return gcd(y, x - y);
}

它的复杂度是 $\mathcal{O}(\max(a,b))$，主要用于高精度运算中，因为模运算（除法）不好实现。

但比如这样：求 $\gcd(1,10^{10000})$ 用更相减损会很慢。除非数据有特殊性质，否则还是老老实实实现高精模运算。

求解 LCM

根据定义，要这么写：

inline int lcm(int x, int y)
{
    return x / gcd(x, y) * y;
}

注意顺序，因为 $x\times y$ 可能会溢出！注意这样的细节，毕竟算法竞赛不是数学竞赛。

质数（素数）

如果一个数 $x(x\in\mathbb{N})$ 的约数仅有 $1$ 和它本身，那么就称 $x$ 是质数（素数），特别地，$0$ 和 $1$ 不是质数，如果一个自然数不是质数，他就是合数。

质数的判定。试除法是最常用的判定质数方法，代码如下：

inline bool is_prime(int x)
{
    if (x < 2) return false;
    for (int i = 2, t = min(x, sqrt(x) + 2); i < t; ++i)
        if (x % i == 0) return false;
    return true; 
}

用反证法可以证明这样做的正确性（即证明合数一定存在一个不超过 $\sqrt{x}$ 的约数）。

“Miller-Rabin”是一个进阶的素数随机判定方法，虽然是随机算法，但没有任何已知数字通过了高级概率性测试被判定为伪素数。感兴趣的读者可以自行了解。

为方便，接下来我们用 $p$ 代表质数。

计算技巧

这是一些编程中的计算技巧。

快速幂

快速幂是一种分治法的应用，但这货在数论中有大量应用。这里直接给出模板。

inline i64 mpow(i64 a, i64 p, i64 m)
{
    a %= p; // 如果 a 在 int 范围内，此步不必要
    i64 res = 1;
    while (p)
    {
        if (p & 1) res = res * a % m;
        a = a * a % m;
        p >>= 1;
    }
    return res % m; // 小心 m = 1
}

快速乘

如果我们要求 64 位整数乘 64 位整数，显然 long long 是要炸掉的。那么用什么？__int128？你觉得毒瘤 CCF 不会使用 128 位整数乘法吗（不用怕，不会的，使用 __int128 就行了，否则接下来的这种方法也不起作用了）。

一种解决方式是把快速幂中的乘法改成加法。但这样实在是太慢了，这里给出一种利用语言特性进行运算的方法，原理请读者自行学习，代码如下：

using i64 = long long;
inline i64 mul(i64 a, i64 b, i64 p)
{
    a %= p, b %= p;
    i64 c = (long double)a * b / p;
    i64 ans = a * b - c * p;
    if (ans < 0) ans += p;
    else if (ans >= p) ans -= p;
    return ans;
}

考场还是能用 __int128 就用。

唯一分解定理

对于一个数 $n$，它可以唯一分解成 $n=\prod p_i^{k_i}$。

这玩意可以证明，我也不证了（这货还叫作算术基本定理，我怎么会证

这玩意的过程也叫做分解质因数。代码如下：

int p[105], c[105], m;
void divide(int n)
{
    m = 0;
    int t = min(n, sqrt(n) + 1);
    for (int i = 2; i <= t; ++i)
        if (n % i == 0)
        {
            p[++m] = i, c[m] = 0;
            while (n % i == 0) n /= i, ++c[m];
        }
    if (n > 1) p[++m] = n, c[m] = 1; // 剩下的这个 n 可能是质数
}

factor 998244353

998244353: 998244353

这货既然敢叫算术基本定理，那么必有过人之处。这里不加证明的给出算术基本定理的常用推论（证明参考《初等数论》）。

给定 $a=\prod_{i=1}^{s} p_i^{\alpha_i}$，

1. $d$ 是 $a$ 的约数的充要条件是 $d=\prod_{i=1}^s p_i^{e_i},e_i\le\alpha_i$，
2. $b=\prod_{i=1}^{s} p_i^{\beta_i}$，那么 $(a,b)=\prod_{i=1}^{s} p_i^{\min(\alpha_i,\beta_i)},[a,b]=\prod_{i=1}^{s} p_i^{\max(\alpha_i,\beta_i)}$，如果我们承认它成立，可以发现 $[a,b](a,b)=ab$ 的原因是 $\min(\alpha,\beta)+\max(\alpha,\beta) =\alpha+\beta$，直呼厉害，
3. 用除数函数 $\tau(a)$（读作 tau，也可记作 $d(a)$）表示 $a$ 的正约数的个数，则 $\tau(a)=\prod_{i=1}^s(\alpha_i+1)$，这也是约数个数定理。其实这是约数条件的推论，对于每个质因子上的幂次，我们取 $0\sim \alpha_i$ 的任意整数，共 $\alpha_i+1$ 个，由乘法原理直接得出，
4. 用除数和函数 $\sigma(a)$（读作 sigma）表示 $a$ 所有正约数的和，则

$$\sigma(a)=\prod_{i=1}^s\left(\sum_{j=0}^{\alpha_i}(p_i)^j\right)=\prod_{i=1}^s\sigma(p_i^{\alpha_i})$$

数论函数

定义域为整数的函数称为数论函数。特别地，对于数论函数 $f$，如果 $\forall a\perp b, f(ab)=f(a)f(b)$，则称 $f$ 为积性函数；如果 $\forall a,b,f(ab)=f(a)f(b)$，则称 $f$ 为完全积性函数。

对于所有积性函数 $f$，有 $n=\prod p_i^{c_i}$，则 $f(n)=\prod f(p_i^{c_i})$。证明很简单，由于第二个式子中的所有 $p_i^{c_i}$ 都必定互质，多次利用积性函数的定义即可证明。

简单数论题

什么，为什么题目来的这么快。
但是不用担心，这里都是一些基础题目，用介绍的知识和一小点小学数学就可以解决。

好，那就来吧^[2]！

线性 GCD

给定 $n(n \le 10^7)$ 个整数，对于每个数 $A_i$，求出删除它以后剩下的数的最大公约数。

还记得 $(a,b,c)=((a,b),c)$ 吗？这个式子告诉我们，GCD 运算存在某种“结合律”，我们可以轻易地把两坨数的 GCD 拼起来。

我们有：

$$Ans_i=(a_1,\cdots,a_{i-1},a_{i+1},\cdots,a_n)=((a_1,\cdots,a_{i-1}),(a_{i+1},\cdots,a_n))$$

处理出前缀 GCD 和后缀 GCD，便可以快速求解了。

[NOIP2001 普及组] 最大公约数和最小公倍数问题

Portal.

有 $p\mid y_0$，我们可以枚举 $y_0$ 的所有约数，然后利用 $(p,q)[p,q]=pq$ 判断是否满足条件即可。

#include <iostream>
#include <cstdio>
#include <algorithm>

using namespace std;

int gcd(int x, int y) { return y == 0 ? x : gcd(y, x % y); }
int lcm(int x, int y) { return x / gcd(x, y) * y; }

int main(void) {
    int x, y, ans = 0;
    cin >> x >> y;
    for (int p = 1; p * p <= y; ++p) 
        if (y % p == 0) {
            if (gcd(p, y / p * x) == x) ++ans;
            if (y / p != p && gcd(y / p, p * x) == x) ++ans;
        }
    cout << ans << endl;
    return 0;
}

[Luogu P1572] 计算分数

Portal.

gcd 可以用于模拟分数运算。

#include <iostream>
#include <cstdio>
#include <string>

#define i64 long long

using namespace std;

i64 gcd(i64 x, i64 y) { return y == 0 ? x : gcd(y, x % y); }
i64 lcm(i64 x, i64 y) { return x / gcd(x, y) * y; }

struct frac {
    i64 p, q;
    frac(i64 p = 0, i64 q = 0) :
        p(p), q(q) {}
    frac operator+ (const frac &a) const {
        frac res;
        res.q = lcm(q, a.q);
        res.p = res.q / q * p + res.q / a.q * a.p;
        i64 g = gcd(res.p, res.q);
        res.p /= g, res.q /= g;
        return res;
    }
};

int main(void) {
    i64 a, b;
    scanf("%lld/%lld", &a, &b);
    frac ans(a, b);
    while (~scanf("%lld/%lld", &a, &b))
        ans = ans + frac(a, b);
    if (ans.q < 0) ans.p = -ans.p, ans.q = -ans.q;
    if (ans.q == 1) printf("%lld\n", ans.p);
    else printf("%lld/%lld\n", ans.p, ans.q);
    return 0;
}

[AHOI2005] 约数研究

Portal.

换一种思路：枚举约数，然后加起来即可。

#include <iostream>
#include <cstdio>

using namespace std;

int main(void) {
    int n, sum = 0;
    cin >> n;
    for (int i = 1; i <= n; i++)
        sum += n / i;
    cout << sum << endl;
    return 0;
}

[Luogu P2651] 添加括号III

Portal.

不难发现除了 $a_2$，其它数都可以放在分子的位置上。然后直接约分，每次用最大公约数除，最后 $a_2=1$ 就是可以。

#include <iostream>
#include <cstdio>

using namespace std;

int gcd(int x, int y) {
    if (y == 0) return x;
    return gcd(y, x % y);
}

int main(void) {
    int T;
    scanf("%d", &T);
    while (T--) {
        int n, a, b;
        scanf("%d%d%d", &n, &a, &b);
        n -= 2;
        b /= gcd(a, b);
        while (n--) {
            scanf("%d", &a);
            b /= gcd(a, b);
        }
        if (b == 1) puts("Yes");
        else puts("No");
    }
    return 0;
}

[Luogu P2660] zzc 种田

Portal.

类似于一个 gcd 的迭代过程。

#include <iostream>
#include <cstdio>

#define i64 long long

using namespace std;

int main(void) {
    i64 a, b, ans = 0;
    cin >> a >> b;
    while (a && b) {
        if (a > b) swap(a, b); // 令 a < b
        ans += a * (b / a) * 4; // 边长为 a 的正方形，种 b / a 个
        b %= a; // 之后 b %= a
    }
    printf("%lld\n", ans);
    return 0;
}

[NOIP2009 提高组] Hankson 的趣味题

Portal.

由于 $[x,b_0]=b_1$，证明 $x$ 是 $b_1$ 的约数，用试除法枚举，判断是否满足条件即可。代码如下：

#include <iostream>
#include <cstdio>
#include <cmath>

using namespace std;

inline int read(void)
{
    int x = 0, c = getchar();
    while (!isdigit(c)) c = getchar();
    while (isdigit(c)) x = (x<<3) + (x<<1) + (c^48), c = getchar();
    return x;
}

int gcd(int x, int y) { return y == 0 ? x : gcd(y, x % y); }
inline int lcm(int x, int y) { return x / gcd(x, y) * y; }

int main(void)
{
    int T = read();
    while (T--)
    {
        int a0 = read(), a1 = read(), b0 = read(), b1 = read(), ans = 0;
        for (int i = 1; i * i <= b1; ++i)
            if (b1 % i == 0)
            {
                if (gcd(i, a0) == a1 && lcm(i, b0) == b1) ++ans;
                if (i != b1 / i)
                {
                    int j = b1 / i;
                    if (gcd(j, a0) == a1 && lcm(j, b0) == b1) ++ans;
                }
            }
        printf("%d\n", ans);
    }
    return 0;
}

[Luogu P1414] 又是毕业季II

Portal.

我们可以计算每个数的约数，然后放到一个数组里面统计。回答的时候从大到小看当前这个约数的数量够不够大。

#include <iostream>
#include <cstdio>
#include <cmath>

using namespace std;

inline int read(void) {
    int x = 0, c = getchar();
    while (!isdigit(c)) c = getchar();
    while (isdigit(c)) x = (x<<3) + (x<<1) + (c^48), c = getchar();
    return x;
}

int n;
int c[1000005];

int main(void) {
    n = read();
    for (int i = 1; i <= n; ++i) {
        int x = read();
        int t = sqrt(x + 0.5);
        for (int j = 1; j <= t; ++j)
            if (x % j == 0) {
                ++c[j];
                if (j * j != x) ++c[x / j];
            }
    }
    for (int i = 1, j = 1000000; i <= n; ++i) {
        while (c[j] < i) --j;
        printf("%d\n", j);
    }
    return 0;
}

[NOIP2013 提高组] 转圈游戏

Portal.

根据题意使用快速幂计算即可。

#include <iostream>
#include <cstdio>

#define i64 long long

using namespace std;

int poww(int a, int p, int m)
{
    int res = 1;
    while (p)
    {
        if (p & 1) res = (i64)res * a % m;
        a = (i64)a * a % m;
        p >>= 1;
    }
    return res % m;
}

int main(void)
{
    int n, m, k, x;
    scanf("%d%d%d%d", &n, &m, &k, &x);
    printf("%d\n", (x + (i64)m * poww(10, k, n)) % n);
    return 0;
}

[NOIP2009 普及组] 细胞分裂

Portal.

将 $m_1$ 分解质因数，然后每个次数乘上 $m_2$，便相当于 $m_1^{m_2}$。对于每一个 $s$，想要让 $m_1^{m_2} \mid s^x$，需要 $m$ 的每个质因子次数都小于 $s^x$ 的。那么依次检查每个质因子，要求 $ex\ge c$ 即可（$e$ 代表将 $s$ 分解质因数后的幂次，$c$ 是 $m$ 处理过后的幂次）。

#include <iostream>
#include <cstdio>

using namespace std;

int n, m;
int m1, m2;
int p[10005], c[10005];

void init(void)
{
    for (int i = 2; i * i <= m1; ++i)
        if (m1 % i == 0)
        {
            p[++m] = i, c[m] = 0;
            while (m1 % i == 0)
            {
                ++c[m];
                m1 /= i;
            }
        }
    if (m1 > 1) p[++m] = m1, c[m] = 1;
    for (int i = 1; i <= m; ++i)
        c[i] *= m2;
}

int main(void)
{
    scanf("%d%d%d", &n, &m1, &m2);
    init();
    int ans = 2e9, s;
    while (n--)
    {
        scanf("%d", &s);
        int x = 0;
        for (int i = 1; i <= m; ++i)
        {
            if (s % p[i] != 0)
            {
                x = 2e9;
                break;
            }
            else
            {
                int res = 0;
                while (s % p[i] == 0)
                {
                    ++res;
                    s /= p[i];
                }
                x = max(x, c[i] / res + (c[i] % res == 0 ? 0 : 1));
            }
        }
        ans = min(ans, x);
    }
    if (ans != 2e9) printf("%d\n", ans);
    else puts("-1");
    return 0;
}

[UVA11889] Benefit

Portal.

当 $a\nmid c$ 的时候无解。否则要求满足下列条件的最小 $x$：

$$\begin{aligned} &\frac{x}{(a,x)}=\frac{c}{a}=d\\ \Longrightarrow & x=dk(k=(a,x))\\ \Longrightarrow & k=(a,dk)\\ \Longrightarrow & \frac{k}{g}=\left(\frac{a}{g},k\times\frac{d}{g}\right)\\ \Longrightarrow & \frac{k}{g}=\left(\frac{a}{g},k\right) \left(\frac{a}{g}\perp\frac{d}{g}\right) \\ \Longrightarrow & \frac{k}{(\frac{a}{g},k)}=g \end{aligned}$$

然后就可以递归求解了。

#include <iostream>
#include <cstdio>

using namespace std;

int gcd(int x, int y) {
    if (y == 0) return x;
    return gcd(y, x % y);
}

int calc(int a, int d) { // b / gcd(a, b) = d
    int g = gcd(a, d);
    if (g == 1) return d;
    return d * calc(a / g, g);
}

int main(void) {
    int T;
    scanf("%d", &T);
    while (T--) {
        int a, c; scanf("%d%d", &a, &c);
        if (c % a != 0) puts("NO SOLUTION");
        else printf("%d\n", calc(a, c / a));
    }   
    return 0;
}

接下来我们会介绍一些较为有趣的数论算法。

素数筛法

素数筛能快速地求出区间内的质数。

模板。

普通筛

对于所有的数，标记所有它的倍数（二倍及以上）为合数。时间复杂度为调和级数 $\mathcal{O}(n\log n)$。代码如下：

int prime[MAXN];
bool vis[MAXN];
for (int i = 2; i <= n; ++i)
{
    if (vis[i] == 0) prime[++tot] = i;
    for (int j = i * 2; j <= n; j += i)
        vis[j] = 1;
}

这种方式很慢，不要在考场上使用。

Eratosthenes 筛法

Eratosthenes 筛法，简称埃氏筛，是算法竞赛中最常用的筛法。

考虑这样一件事情：如果 $x$ 是合数，那么 $x$ 的倍数也一定是合数，它们肯定会被某个相同的质因子筛掉，那我们就没有必要筛合数了。
如果我们从小到大考虑每个数，然后同时把当前这个数的所有（比自己大的）倍数记为合数，那么运行结束的时候没有被标记的数就是素数了。

代码如下：

int prime[MAXN];
bool vis[MAXN];
for (int i = 2; i <= n; ++i)
    if (vis[i] == 0)
    {
        prime[++tot] = i;
        for (int j = i * 2; j <= n; j += i)
            vis[j] = 1;
    }

埃氏筛的复杂度是 $\mathcal{O}(\sum_{p\leqslant n} \cfrac{n}{p}=n \log\log n)$，效率较高，但仍需较好的 IO 优化才可通过刚才的模板。

这个代码还可以加一些优化：内层循环不必从 $i\times 2$ 开始，它在筛 $2$ 的时候就已经被筛掉了；外层循环也不必到 $n$ 结束，需要筛的已经在 $\sqrt{n}$ 前筛掉了（但这道题中这个两个都不能用，因为分别会爆 int 和无法统计素数）。但实际上，优化后的复杂度并没有改变。

埃氏筛的速度已经足够快，如果不是时间特别吃紧，考场上使用埃氏筛即可。

刚才是“质数的整数倍是合数“，还有一种思路：整数的质数倍是合数。代码如下：

for (int i = 2; i <= n; ++i) {
    if (!vis[i]) prime[++tot] = i;
    for (int j = 1; j <= tot && i * prime[j] <= n; ++j) 
        vis[i * prime[j]] = true;
}

欧拉筛

欧拉筛（也称作线性筛）的时间复杂度为 $\mathcal{O}(n)$，是常用筛法中最快的筛法。

埃氏筛的问题在于它仍然会重复标记合数（从 $x^2$ 开始），我们需要设计一个算法，使得合数只会被它最小的质因子筛一次，这就是欧拉筛（由于欧拉筛的这个特性，它也常被用做求合数最小质因子的方法）。

我们设一个 $v$ 数组代表每个数的最小质因子，设质数集合为 $p$，那么可以这样维护 $v$：

1. 一次考虑 $2\sim n$ 间的每一个数 $i$。
2. 如果 $v[i]$ 没有值，那么说明它是素数，把它加进 $p$ 中。
3. 接下来对于每个质数 $p[j]$，令 $v[i\times p[j]]=p[j]$。由于 $p[j] \le v[i]$，所以 $p[j]$ 一定是 $i\times p[j]$ 的最小质因子。

由于所有的合数只会被标记一次，时间复杂度为 $\mathcal{O}(n)$，严谨的复杂度和正确性证明较为复杂，感兴趣的读者可以自行了解。代码如下：

int n, tot = 0;
int v[MAXN], prime[MAXN];
inline void primes(void)
{
    for (int i = 2; i <= n; ++i)
    {
        if (v[i] == 0) // 是质数
        {
            v[i] = i;
            prime[++tot] = i;
        }
        for (int j = 1; j <= tot && i * prime[j] <= n; ++j) // 给 i 乘上一个质数，乘积小于 n 才循环
        {
            if (prime[j] > v[i]) break; // i 的最小质因子比 prime[j] 更小，break
            v[i * prime[j]] = prime[j]; // 标记
        }
    }
}

实际上如果不用记录最小值因子的话，也可以这么做：

for (int i = 2; i <= n; ++i) {
    if (!vis[i]) prime[++tot] = i;
    for (int j = 1; j <= tot && i * prime[j] <= n; ++j) {
        vis[i * prime[j]] = true;
        if (i % prime[j] == 0) break;
    }
}

只多了一个 break，就让埃氏筛变成了欧拉筛，每个数只被筛了一次。

区间筛

[Luogu P1835] 素数密度.

给定区间 $[L,R](1\le L\le R < 2^{31},R-L\le 10^6)$，求区间中的素数个数。

$L,R$ 这么大，但是 $R-L$ 的值却这么小，而且任何一个合数 $n$ 都一定包含一个不超过 $\sqrt{n}$ 的质因子，那么我们只需要用筛法求出 $2\sim \sqrt{R}$ 以内的所有质数，然后筛掉这些数的倍数即可。如果采用线性筛，那么时间复杂度为 $O(\sqrt{r}+(r-l+1)\log\log \sqrt{r})$。代码如下：

#include <iostream>
#include <cstdio>
#include <cmath>
#define LL long long

using namespace std;

LL L, R;
bool is_prime[50000];
bool flag[1000005];

inline void solve(int n)
{
    for (int i = 2; i <= n; ++i) is_prime[i] = 1;
    int g = sqrt(n) + 1;
    for (int i = 2; i <= g; ++i) if (is_prime[i])
        for (int j = i * i; j <= n; j += i)
            is_prime[j] = 0;
}

int main(void)
{
    cin >> L >> R;
    if (L == 1)
        L = 2;
    if (L > R)
    {
        puts("0");
        return 0;
    }
    int m = sqrt(R) + 1;
    solve(m);

    for (int i = 0; i <= int(R - L); ++i)
        flag[i] = 1;
    for (int i = 2; i <= m; ++i)
        if (is_prime[i])
        {
            LL j = L;
            while (j % i != 0)
                ++j;
            for (; j <= R; j += i)
                if (i != j)
                    flag[j - L] = 0;
        }
    int ans = 0;
    for (int i = 0; i <= R - L; ++i)
        if (flag[i])
            ++ans;
    cout << ans << endl;
    return 0;
}

筛法与积性函数

筛法的作用远不止筛出素数，它们还可以在筛素数的同时求出积性函数。在欧拉函数一节中，我们将看到如何用筛法求解欧拉函数；在《数论进阶》中，我们还会更深入地探讨筛法的应用。

扩展欧几里得算法

扩展欧几里得算法用于求解不定方程 $ax+by=\gcd(a,b)$ 的整数解，是数论基础算法之一。

裴蜀定理

裴蜀定理。$\forall a,b, ~\exists x,y,~\text{i.e.}~ ax+by=\gcd(a,b)$

关于裴蜀定理的证明，网上资料较多，这里不做赘述。

模板。

根据裴蜀定理，$ax+by=\gcd(a,b)$，也就是说 $ax+by>0$ 时，最小为 $\gcd(a,b)$，也就是说答案是 $\gcd(A)$，而且是一个正数。

#include <iostream>
#include <cstdio>

using namespace std;

int gcd(int x, int y)
{
    if (y == 0) return x;
    return gcd(y, x % y);
}

int main(void)
{
    int x, ans = 0;
    scanf("%d", &x);
    while (scanf("%d", &x) == 1) ans = gcd(ans, x < 0 ? -x : x);
    printf("%d\n", ans);
    return 0;
}

扩展欧几里得

由于 $ax+by=\gcd(a,b)=\gcd(b,a\bmod b)$ 且 $ax+by=(\lfloor\cfrac{a}{b}\rfloor\times b+a\bmod b)x+by=b(y+\lfloor\cfrac{a}{b}\rfloor x)+(a\bmod b)x$，所以求不定方程 $ax+by=\gcd(a,b)$，相当于求 $b(y+\lfloor\cfrac{a}{b}\rfloor x)+(a\bmod b)x=\gcd(b, a\bmod b)$。

当然这里把代码实现考虑了进去（还记得为什么求 gcd 要把 $x$ 和 $y$ 写反吗），所以新式子中的 $(y+\lfloor\cfrac{a}{b}\rfloor x)$ 相当于 $y'$，$x$ 相当于 $x'$。

进而有 $y = y' - \lfloor\cfrac{a}{b}\rfloor\times x'$，$x$ 是 $x'$。

所以：

void exgcd(int a, int b, int &x, int &y) {
    if (b == 0) return x = 1, y = 0, void(); // 显然
    exgcd(b, a%b, y, x); // 根据刚才推出的结论
    // 此时已经求出了 x', y'
    y -= a / b * x;
}

那么如何求出方程的所有解呢？

设用 exgcd 求出的解为 $x_0,y_0$，那么显然有 $x=x_0+k\times\cfrac{b}{\gcd(a,b)},y=y_0-k\times\cfrac{a}{\gcd(a,b)}$。

最小正数 $x$ 显然为 $x_0 \bmod \cfrac{b}{\gcd(a,b)}$。

模板。

根据裴蜀定理，$ax+by$ 的答案只能是 $\gcd(a,b)$ 的倍数，也就是需要满足 $(a,b)\mid c$。直接使用 exgcd 算法，求出来的解，然后再乘上 $c\div (a,b)$，就能够得到一组原方程的解。

具体实现见代码：

// 笔记原因，只摘录了重要代码
i64 exgcd(i64 a, i64 b, i64 &x, i64 &y) {
    if (b == 0) return x = 1, y = 0, a;
    i64 d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

int main(void) {
    i64 a = read(), b = read(), c = read(), x, y;
    i64 g = exgcd(a, b, x, y);
    if (c % g != 0) puts("-1");
    else {
        x *= c / g, y *= c / g;
        i64 p = b / g, q = a / g, k;
        if (x < 0) k = ceil((1.0 - x) / p), x += p * k, y -= q * k;
        else if (x >= 0) k = (x - 1) / p, x -= p * k, y += q * k;
        if (y > 0) {
            printf("%lld %lld %lld ", (y - 1) / q + 1, x, (y - 1) % q + 1);
            printf("%lld %lld\n", x + (y - 1) / q * p, y)
        }
        else printf("%lld %lld\n", x, y + q * (i64)ceil((1.0 - y) / q));
    }
}

欧拉函数

欧拉函数，即 $\varphi(n)$，代表 $1\sim n$ 中与 $n$ 互质的数的个数。

比如在 $1$~$12$ 中，$1,5,7,11$ 这四个数与 $12$ 互质，所以 $\varphi(12)=4$。

求解

怎么求欧拉函数呢？设 $n=\prod p_i^{k_i}$，则有：

$$\varphi(n)=n\prod(1-\cfrac{1}{p_i})=n\prod \cfrac{p_i-1}{p_i}$$

证明：

根据结论有如下代码：

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

时间复杂度为 $\mathcal{O}(\sqrt{n})$，当然你可以用 Pollard Rho 来加速，不过没必要。

性质

1. $\forall n > 1, 1\sim n$ 中与 $n$ 互质的数的和为 $n\times \varphi(n) \div 2$。且 $2\mid\varphi(n)(n>2)$。