ChungZH修改

2023-07-29 08:00:00

最大流

Luogu-P1231 教辅的组成

三倍经验： Luogu-P1402 酒店之王 / Luogu-P2891 [USACO07OPEN] Dining G

一眼丁真建图：S->练习册->书->答案->T

然而是错的。很明显，书有可能被多次匹配，与题意不符。

正确的建图：S->练习册->书（拆点）->答案->T

为什么中间层的书要拆点呢？因为一本书不能被重复选用。我们的目的是保证一本书流出的流量只能是 $1$。所以我们把每个代表书的点拆成两个点，左边的点和练习册连边，右边的点和答案连边；左右对应点之间也要连一条容量为 $1$ 的边。

Luogu-P2764 最小路径覆盖问题

定理：最小路径覆盖数=$|G|$-二分图最大匹配数

首先我们假设现在原图内每个点都是一条路径，此时最少路径数为 $n$。

考虑合并路径，当且仅当两条路径首尾相连的时候可以合并。

将点 $x$ 拆成出点 $x$ 和入点 $x+n$，当我们连接 $u, v$ 时，转化为连接 $u, v+n$。将 $S$ 与所有 $u$ 连边，将所有 $u+n$ 与 $T$ 连边。所有边的容量都为 $1$。

在一开始每个点都是一条独立的路径，每次合并将两条路径合并为一条路径，那么最终路径即为点数减去最大匹配数，这样求得的路径覆盖即为最小路径覆盖。

对于输出路径，用 to 记录下一个节点，tag 标记该节点前面是否还有点。

 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
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

ll dfs(int u, int t, ll flow) {
 if (u == t) return flow;
 ll ans = 0;
 vis[u] = 1;
 for (int &i = cur[u]; ~i; i = e[i].nxt) { // i 用引用：当前弧优化
 int v = e[i].to;
 if (!vis[v] && dep[v] == dep[u] + 1 && e[i].cap > e[i].flow) {
 ll d = dfs(v, t, min(flow - ans, e[i].cap - e[i].flow));
 if (!d) continue;
 ans += d;
 e[i].flow += d;
 e[i ^ 1].flow -= d;
 to[u] = v;
 if (u != S) { tag[v - n] = 1; }
 if (ans == flow) break; // 剪枝，残余流量用尽，停止增广
 }
 }
 vis[u] = 0;
 return ans;
}
ll Dinic(int s, int t) {
 ll ret = 0;
 while (bfs(s, t)) {
 memcpy(cur, fir, sizeof cur);
 ret += dfs(s, t, INF);
 }
 for (int i = 1; i <= n; i++) {
 if (tag[i] == 0) {
 cout << i << " ";
 int x = i;
 while (to[x] && to[x] != t) {
 cout << to[x] - n << " ";
 x = to[x] - n;
 }
 cout << "\n";
 }
 }
 return ret;
}
int main() {
 ios::sync_with_stdio(false);

 memset(fir, -1, sizeof fir);
 cin >> n >> m;
 S = 2 * n + 1, T = S + 1;
 for (int i = 0; i < m; i++) {
 int x, y;
 cin >> x >> y;
 addEdge(x, y + n, 1);
 }
 for (int i = 1; i <= n; i++) {
 addEdge(S, i, 1);
 addEdge(i + n, T, 1);
 }
 int t = Dinic(S, T);
 cout << n - t << endl;
 return 0;
}

最小割

方格取数问题

对矩阵进行黑白染色，那么比如选了一个黑格，那么与它相邻的白格不能选。可以计算不选的方格，其余的都选。容易想到最小割。

那么建图就很容易了。对于每个点 $(i,j)$，数值为 $a_{i,j}$，如果是白色点，那么从源点向它连一条容量 $a_{i,j}$ 的边；否则就从这个点向汇点连一条容量 $a_{i,j}$ 的边。这样我们就把不选的点转化成了割掉的边。为了保证相邻点直接的边不会被割掉，我们从每个白点向与其相邻的黑点连容量为 $INF$ 的边（他们之间的边只从白连到黑，避免重复计算）。答案为全局和减去最小割。

Luogu-P2598 [ZJOI2009] 狼和羊的故事

原点向所有狼连容量 $INF$ 的边，所有羊向汇点连容量 $INF$ 的边，所有点向四周连容量为 $1$ 的边。

最小割即为答案，因为狼和羊之间的边都被割掉了。

Luogu-P1361 小M的作物

二者取其一模型。

考虑点集 ${u, v, w}$ 对集合 $A$ 的贡献，加入其中一者被割进 $B$，那这个点集就没有贡献，即点集与 $A$ 的连边要割掉。

用一个虚点 $x$ 从 $S$ 连一条边代表贡献，如果其中一个点被割进了集合 $B$，那这条代表贡献的边也要割掉，而 $x$ 到 $u, v, w$ 的边不能断开。所以 $(S, x)$ 的容量为 $c$，$(x, u), (x, v), (x, w)$ 的容量都是 $INF$（保证不被断开）。

答案为总收益减去最小割。

Luogu-P2762 太空飞行计划问题

两倍经验：Luogu-P3410 拍照

最大权值闭合图模型。

Luogu-P4177 [CEOI2008] order

如果没有租机器，那这道题就是纯最大权值闭合图模型。

图中对于工作和它所需的机器之间边的容量为 $INF$，所以这条边不可能被割，意义即为选择这个工作就必须购买这个机器。那么租用就可以将这条边的容量改为 $b_{ij}$，可以用 $b_{ij}$ 割掉这条边，表示选择这个工作后，可以花费 $b_{ij}$ 的代价代替购买机器。

2023-07-29 08:00:00

最大流

概念

• 容量：$capacity(e)$ 表示一条有向边 $e(u, v)$ 的最大允许的流量。
• 流量：$flow(e)$ 表示一条有向边 $e(u, v)$ 总容量中已被占用的流量。
• 剩余流量（残量）：$w(e) = c(e)-f(e)$，表示当前时刻某条有向边 $e(u, v)$ 总流量中未被占用的部分。
• 反向边：原图中每一条有向边在残量网络中都有对应的反向边，反向边的容量为 $0$，容量的变化与原边相反；『反向边』的概念是相对的，即一条边的反向边的反向边是它本身。
• 残量网络：在原图的基础之上，添加每条边对应的反向边，并储存每条边的当前流量。残量网络会在算法进行的过程中被修改。
• 增广路（augmenting path）：残量网络中从源点到汇点的一条路径，增广路上所有边中最小的剩余容量为增广流量。
• 增广（augmenting）：在残量网络中寻找一条增广路，并将增广路上所有边的流量加上增广流量的过程。
• 层次：$level(u)$ 表示节点 $u$ 在层次图中与源点的距离。
• 层次图：在原残量网络中按照每个节点的层次来分层，只保留相邻两层的节点的图，满载（即流量等于容量）的边不存在于层次图中。

流的三个重要性质：

1. 容量限制：对于每条边，流经该边的流量不得超过该边的容量，即，$f(u,v)\leq c(u,v)$
2. 斜对称性：每条边的流量与其相反边的流量之和为 0，即 $f(u,v)=-f(v,u)$
3. 流守恒性：从源点流出的流量等于汇点流入的流量，即 $\forall x\in V-{s,t},\sum_{(u,x)\in E}f(u,x)=\sum_{(x,v)\in E}f(x,v)$

最大流问题：指定合适的流 $f$，以最大化整个网络的流量（即 $\sum_{(s,v)\in E}f(s,v)$）。

Ford-Fulkerson 增广

增广路指一条从 $s$ 到 $t$ 的路径，路径上每条边残余容量都为正。把残余容量为正（$w(u, v) \gt 0$）的边设为可行边，然后搜索寻找增广路。

找到一条增广路后，这条路能够增广的流值由路径上边的最小残留容量 $w(u, v)$（记为 $flow$）决定。将这条路径上每条边的 $w(u, v)$ 减去 $flow$。最后，路径上每条边的反向边残留容量要加上 $flow$。

为什么增广路径上每条边的反向边的残留容量值要加上 $flow$？因为斜对称性，残量网络=容量网络-流量网络，容量网络不变时，流量网络上的边的流量增加 $flow$，反向边流量减去 $flow$，残量网络就会发生相反的改变。从另一个角度来说，这个操作可以理解为「退流」，给了我们一个反悔的机会，让增广路的顺序不受限制。

增广路算法好比是自来水公司不断的往水管网里一条一条的通水。

这个算法基于増广路定理：网络达到最大流当且仅当残留网络中没有増广路。

建图技巧：从 $2$ 开始对边编号，这样 $i$ 的反向边就是 $i \oplus 1$。

Edmonds–Karp 算法

用 BFS 找增广路：

• 如果在 $G_f$ 上我们可以从 $s$ 出发 BFS 到 $t$，则我们找到了新的增广路。
• 对于增广路 $p$，我们计算出 $p$ 经过的边的剩余容量的最小值 $\Delta = \min_{(u, v) \in p} c_f(u, v)$。我们给 $p$ 上的每条边都加上 $\Delta$ 流量，并给它们的反向边都退掉 $\Delta$ 流量，令最大流增加了 $\Delta$。
• 因为我们修改了流量，所以我们得到新的 $G_f$，我们在新的 $G_f$ 上重复上述过程，直至增广路不存在，则流量不再增加。

显然，单轮 BFS 增广的时间复杂度是 $O(m)$。增广总轮数的上界是 $O(nm)$。相乘后得到 EK 算法的时间复杂度是 $O(nm^2)$

代码实现中，$w$ 表示边的容量。

 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
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int MAXN = 250, MAXM = 5005;
const ll INF = 0x3f3f3f3f;

struct Edge {
 int from, to;
 ll w;
} edges[2 * MAXM];
int n, m, s, t, tot = 1;
vector<int> G[MAXN]; // G: x 的所有边在 edges 中的下标
ll a[MAXN]; // a: BFS 过程中最近接近点 x 的边给它的最大流
int p[MAXN]; // p: 最近接近 x 的边
void addEdge(int from, int to, ll w) {
 edges[++tot] = {from, to, w};
 edges[++tot] = {to, from, 0};
 G[from].push_back(tot - 1);
 G[to].push_back(tot);
}
bool bfs(int s, int t) {
 memset(a, 0, sizeof a);
 memset(p, 0, sizeof p);
 queue<int> Q;
 Q.push(s);
 a[s] = INF;
 while (!Q.empty()) {
 int x = Q.front();
 Q.pop();
 for (int i = 0; i < G[x].size(); i++) {
 Edge &e = edges[G[x][i]];
 if (!a[e.to] && e.w) {
 p[e.to] = G[x][i]; // G[x][i] 是最接近 e.to 的边
 a[e.to] = min(a[x], e.w); // 最接近 e.to 的边赋给它的流
 Q.push(e.to);
 if (e.to == t) return true;
 }
 }
 }
 return false;
}
ll EK(int s, int t) {
 ll flow = 0;
 while (bfs(s, t)) {
 for (int u = t; u != s; u = edges[p[u]].from) { // 回溯路径
 edges[p[u]].w -= a[t];
 edges[p[u] ^ 1].w += a[t];
 }
 flow += a[t];
 }
 return flow;
}
int main() {
 ios::sync_with_stdio(false);
 cin >> n >> m >> s >> t;
 for (int i = 0; i < m; i++) {
 int u, v;
 ll w;
 cin >> u >> v >> w;
 addEdge(u, v, w);
 }
 cout << EK(s, t) << endl;
 return 0;
}

Dinic 算法

EK 在有些情况中比较低效，我们引入另一个算法：Dinic。

考虑在增广前先对 $G_f$ 做 BFS 分层，即根据结点 $u$ 到源点 $s$ 的距离 $d(u)$ 把结点分成若干层。令经过 $u$ 的流量只能流向下一层的结点 $v$，即删除 $u$ 向层数标号相等或更小的结点的出边，我们称 $G_f$ 剩下的部分为层次图（Level Graph）。形式化地，我们称 $G_L = (V, E_L)$ 是 $G_f = (V, E_f)$ 的层次图，其中 $E_L = \left{ (u, v) \mid (u, v) \in E_f, d(u) + 1 = d(v) \right}$。

如果我们在层次图 $G_L$ 上找到一个最大的增广流 $f_b$，使得仅在 $G_L$ 上是不可能找出更大的增广流的，则我们称 $f_b$ 是 $G_L$ 的阻塞流（Blocking Flow）。

定义层次图和阻塞流后，Dinic 算法的流程如下。

1. 在 $G_f$ 上 BFS 出层次图 $G_L$。
2. 在 $G_L$ 上 DFS 出阻塞流 $f_b$。
3. 将 $f_b$ 并到原先的流 $f$ 中，即 $f \leftarrow f + f_b$。
4. 重复以上过程直到不存在从 $s$ 到 $t$ 的路径。

此时的 $f$ 即为最大流。

我们还需要当前弧优化。注意到在 $G_L$ 上 DFS 的过程中，如果结点 $u$ 同时具有大量入边和出边，并且 $u$ 每次接受来自入边的流量时都遍历出边表来决定将流量传递给哪条出边，则 $u$ 这个局部的时间复杂度最坏可达 $O(|E|^2)$。事实上，如果我们已经知道边 $(u, v)$ 已经增广到极限（边 $(u, v)$ 已无剩余容量或 $v$ 的后侧已增广至阻塞），则 $u$ 的流量没有必要再尝试流向出边 $(u, v)$。据此，对于每个结点 $u$，我们维护 $u$ 的出边表中第一条还有必要尝试的出边。习惯上，我们称维护的这个指针为当前弧，称这个做法为当前弧优化。

将单轮增广的时间复杂度 $O(nm)$ 与增广轮数 $O(n)$ 相乘，Dinic 算法的时间复杂度是 $O(n^2m)$。

对于 Dinic 算法的复杂度，有如下 $3$ 种情况：

• 一般的网络图：$O(n^2m)$
• 单位容量的图：$O(\min(\sqrt m,n^{\frac{2}{3}})\cdot m)$
• 二分图：$O(m\sqrt n)$

（Dinic 中 DFS 最好要写 vis，因为求费用流一定要）

错误：有一次漏判了 dep[v] == dep[u] + 1！玄学：若 dep 默认值为 $0$，则一定要 $dep[S]=1$！！！

 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
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77

#include <bits/stdc++.h>
using namespace std;

using ll = long long;
const int MAXN = 250, MAXM = 5005;
const ll INF = 0x3f3f3f3f;

struct Edge {
 int to, nxt;
 ll cap, flow;
} e[2 * MAXM];
int fir[MAXN];
int n, m, s, t, tot = 1;
int dep[MAXN], cur[MAXN]; // cur 记录当前弧
bool vis[MAXN];
void addEdge(int from, int to, ll w) {
 e[++tot] = {to, fir[from], w, 0};
 fir[from] = tot;
 e[++tot] = {from, fir[to], 0, 0};
 fir[to] = tot;
}
bool bfs(int s, int t) {
 queue<int> Q;
 memset(dep, 0, sizeof dep);
 dep[s] = 1;
 Q.push(s);
 while (!Q.empty()) {
 int x = Q.front();
 Q.pop();
 for (int i = fir[x]; ~i; i = e[i].nxt) {
 int v = e[i].to;
 if (!dep[v] && e[i].cap > e[i].flow) {
 dep[v] = dep[x] + 1;
 Q.push(v);
 }
 }
 }
 return dep[t];
}
ll dfs(int u, int t, ll flow) {
 if (u == t) return flow;
 ll ans = 0;
 vis[u] = 1;
 for (int &i = cur[u]; ~i; i = e[i].nxt) { // i 用引用：当前弧优化
 int v = e[i].to;
 if (!vis[v] && dep[v] == dep[u] + 1 && e[i].cap > e[i].flow) {
 ll d = dfs(v, t, min(flow - ans, e[i].cap - e[i].flow));
 ans += d;
 e[i].flow += d;
 e[i ^ 1].flow -= d;
 if (ans == flow) break; // 剪枝，残余流量用尽，停止增广
 }
 }
 vis[u] = 0;
 return ans;
}
ll Dinic(int s, int t) {
 ll ans = 0;
 while (bfs(s, t)) {
 memcpy(cur, fir, sizeof cur); // 当前弧优化
 ans += dfs(s, t, INF);
 }
 return ans;
}
int main() {
 ios::sync_with_stdio(false);
 cin >> n >> m >> s >> t;
 memset(fir, -1, sizeof fir);
 for (int i = 0; i < m; i++) {
 int u, v;
 ll w;
 cin >> u >> v >> w;
 addEdge(u, v, w);
 }
 cout << Dinic(s, t) << endl;
 return 0;
}

最小割

概念

• 割：对于一个网络流图 $G=(V,E)$，其割的定义为一种 点的划分方式：将所有的点划分为 $S$ 和 $T=V-S$ 两个集合，其中源点 $s\in S$，汇点 $t\in T$。
• 割的容量：我们的定义割 $(S,T)$ 的容量 $c(S,T)$ 表示所有从 $S$ 到 $T$ 的边的容量之和，即 $c(S,T)=\sum_{u\in S,v\in T}c(u,v)$。当然我们也可以用 $c(s,t)$ 表示 $c(S,T)$。
• 最小割：最小割就是求得一个割 $(S,T)$ 使得割的容量 $c(S,T)$ 最小。

最大流最小割定理

证明：

定理：$f(s,t){\max}=c(s,t){\min}$

对于任意一个可行流 $f(s,t)$ 的割 $(S,T)$，我们可以得到：

$$ f(s,t)=S\text{出边的总流量}-S\text{入边的总流量}\le S\text{出边的总流量}=c(s,t) $$

如果我们求出了最大流 $f$，那么残余网络中一定不存在 $s$ 到 $t$ 的增广路经，也就是 $S$ 的出边一定是满流，$S$ 的入边一定是零流，于是有：

$$ f(s,t)=S\text{出边的总流量}-S\text{入边的总流量}=S\text{出边的总流量}=c(s,t) $$

结合前面的不等式，我们可以知道此时 $f$ 已经达到最大。

实现

最小割：最大流

方案：我们可以通过从源点 $s$ 开始 DFS，每次走残量大于 $0$ 的边，找到所有 $S$ 点集内的点。

最小割边数量：如果需要在最小割的前提下最小化割边数量，那么先求出最小割，把没有满流的边容量改成 $\infty$，满流的边容量改成 $1$，重新跑一遍最小割就可求出最小割边数量；如果没有最小割的前提，直接把所有边的容量设成 $1$，求一遍最小割就好了。

模型一：二者取其一问题

有 $n$ 个物品和两个集合 $A,B$，如果一个物品没有放入 $A$ 集合会花费 $a_i$，没有放入 $B$ 集合会花费 $b_i$；还有若干个形如 $u_i,v_i,w_i$ 限制条件，表示如果 $u_i$ 和 $v_i$ 同时不在一个集合会花费 $w_i$。每个物品必须且只能属于一个集合，求最小的代价。

这是一个经典的 二者选其一 的最小割题目。我们对于每个集合设置源点 $s$ 和汇点 $t$，第 $i$ 个点由 $s$ 连一条容量为 $a_i$ 的边、向 $t$ 连一条容量为 $b_i$ 的边。对于限制条件 $u,v,w$，我们在 $u,v$ 之间连容量为 $w$ 的双向边。

注意到当源点和汇点不相连时，代表这些点都选择了其中一个集合。如果将连向 $s$ 或 $t$ 的边割开，表示不放在 $A$ 或 $B$ 集合，如果把物品之间的边割开，表示这两个物品不放在同一个集合。

最小割就是最小花费。

模型二：最大权值闭合图

最大权值闭合图，即给定一张有向图，每个点都有一个权值（可以为正或负或 $0$），你需要选择一个权值和最大的子图，使得子图中每个点的后继都在子图中。

做法：建立超级源点 $s$ 和超级汇点 $t$，若节点 $u$ 权值为正，则 $s$ 向 $u$ 连一条有向边，边权即为该点点权；若节点 $u$ 权值为负，则由 $u$ 向 $t$ 连一条有向边，边权即为该点点权的相反数。原图上所有边权改为 $\infty$。跑网络最大流，将所有正权值之和减去最大流，即为答案。

几个小结论来证明：

1. 每一个符合条件的子图都对应流量网络中的一个割。因为每一个割将网络分为两部分，与 $s$ 相连的那部分满足没有边指向另一部分，于是满足上述条件。这个命题是充要的。
2. 最小割所去除的边必须与 $s$ 和 $t$ 其中一者相连。因为否则边权是 $\infty$，不可能成为最小割。
3. 我们所选择的那部分子图，权值和 $=$ 所有正权值之和 $-$ 我们未选择的正权值点的权值之和 $+$ 我们选择的负权值点的权值之和。当我们不选择一个正权值点时，其与 $s$ 的连边会被断开；当我们选择一个负权值点时，其与 $t$ 的连边会被断开。断开的边的边权之和即为割的容量。于是上述式子转化为：权值和 $=$ 所有正权值之和 $-$ 割的容量。
4. 于是得出结论，最大权值和 $=$ 所有正权值之和 $-$ 最小割 $=$ 所有正权值之和 $-$ 最大流。

P2057 [SHOI2007] 善意的投票 / [JLOI2010] 冠军调查 Title

小M的作物Title

最小费用最大流

费用流

给定一个网络 $G=(V,E)$，每条边除了有容量限制 $c(u,v)$，还有一个单位流量的费用 $w(u,v)$。

当 $(u,v)$ 的流量为 $f(u,v)$ 时，需要花费 $f(u,v)\times w(u,v)$ 的费用。

$w$ 也满足斜对称性，即 $w(u,v)=-w(v,u)$。

则该网络中总花费最小的最大流称为 最小费用最大流，即在最大化 $\sum_{(s,v)\in E}f(s,v)$ 的前提下最小化 $\sum_{(u,v)\in E}f(u,v)\times w(u,v)$。

SSP 算法

SSP（Successive Shortest Path）算法是一个贪心的算法。它的思路是每次寻找单位费用最小的增广路进行增广，直到图上不存在增广路为止。

如果图上存在单位费用为负的圈，SSP 算法无法正确求出该网络的最小费用最大流。此时需要先使用消圈算法消去图上的负圈。

时间复杂度： 如果使用 Bellman–Ford 求解最短路，每次找增广路的时间复杂度为 $O(nm)$。设该网络的最大流为 $f$，则最坏时间复杂度为 $O(nmf)$。事实上，SSP 算法是伪多项式时间的。

实现

可以在 Dinic 基础上改进，将 BFS 求分层图换为 SPFA（有负权边不能用 Dijkstra）求一条单位费用之和最小的路径，也就是把 $w(u, v)$ 当做边权然后在残量网络上求最短路，在 DFS 中一定要用 vis 数组记录点是否访问过。这样就可以求得最小费用最大流了。

建反向边时，费用取反即可。

 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
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

#include <bits/stdc++.h>
using namespace std;

using ll = long long;
const int MAXN = 100055, MAXM = 100005;
const int INF = 0x3f3f3f3f;

struct Edge {
 int to, nxt, cap, cost;
} e[2 * MAXM];
int fir[MAXN];
int n, m, s, t, tot = 1;
int dis[MAXN];
bool vis[MAXN];
int cur[MAXN]; // cur 记录当前弧
void addEdge(int from, int to, int w, int cost) {
 e[++tot] = {to, fir[from], w, cost};
 fir[from] = tot;
 e[++tot] = {from, fir[to], 0, -cost};
 fir[to] = tot;
}
bool spfa(int s, int t) {
 memset(dis, 0x3f, sizeof dis);
 queue<int> q;
 q.push(s), dis[s] = 0, vis[s] = 1;
 while (!q.empty()) {
 int f = q.front();
 q.pop();
 vis[f] = 0;
 for (int i = fir[f]; ~i; i = e[i].nxt) {
 int v = e[i].to;
 if (e[i].cap && dis[f] + e[i].cost < dis[v]) {
 dis[v] = dis[f] + e[i].cost;
 if (!vis[v]) {
 vis[v] = 1;
 q.push(v);
 }
 }
 }
 }
 return dis[t] < INF;
}
int dfs(int u, int t, int flow) {
 if (u == t) return flow;
 int ans = 0;
 vis[u] = 1;
 for (int &i = cur[u]; ~i; i = e[i].nxt) { // i 用引用：当前弧优化
 int v = e[i].to;
 if (!vis[v] && dis[v] == dis[u] + e[i].cost && e[i].cap) {
 int d = dfs(v, t, min(flow - ans, e[i].cap));
 ans += d;
 e[i].cap -= d;
 e[i ^ 1].cap += d;
 // 易错：这里不写 return ans，应为 break，需要执行下面的 vis[u] = 0
 if (ans == flow) break; // 剪枝，残余流量用尽，停止增广
 }
 }
 vis[u] = 0;
 return ans;
}
pair<int, int> Dinic(int s, int t) {
 int maxFlow = 0, minCost = 0;
 while (spfa(s, t)) {
 memcpy(cur, fir, sizeof cur); // 当前弧优化
 int x = dfs(s, t, INF);
 maxFlow += x;
 minCost += x * dis[t];
 }
 return make_pair(maxFlow, minCost);
}
int main() {
 cin >> n >> m >> s >> t;
 memset(fir, -1, sizeof fir);
 for (int i = 0; i < m; i++) {
 int u, v, w, c;
 cin >> u >> v >> w >> c;
 addEdge(u, v, w, c);
 }
 pair<int, int> ans = Dinic(s, t);
 cout << ans.first << " " << ans.second << endl;
 return 0;
}

参考资料

• OI Wiki
• 最大权闭合子图 - kikokiko
• 网络流 - Siyuan
• 最详细（也可能现在不是了）网络流建模基础 - ~victorique~

2023-05-21 08:00:00

CF-559C Gerald and Giant Chess / AtCoder DP-Y Grid 2

给定一个 $H*W$ 的棋盘，棋盘上只有 $N$ 个格子是黑色的，其他格子都是白色的。在棋盘左上角有一个卒，每一步可以向右或者向下移动一格，并且不能移动到黑色格子中。求这个卒从左上角移动到右下角，一共有多少种可能的路线。

$(1 ≤ h, w ≤ 105, 1 ≤ n ≤ 2000)$

$O(hw)$ 的暴力 DP 很好想，但是过不了。

假设没有障碍，从 $(1, 1)$ 到 $(i, j)$ 的方案数是 $C_{i+j-2}^{i-1}$（等于 $C_{i+j-2}^{j-1}$）。可以这么理解：可以用 $D, R$ 来表示一条路径，那么从 $(1, 1)$ 到 $(i, j)$ 的路径中有 $i-1$ 个 $D$ 和 $j-1$ 个 $R$。于是问题转化为从 $i+j-2$ 个位置中选 $i-1$ 个放 $D$ 的方案数。

如果有一个障碍，从正面统计方案数很困难，正难则反，考虑将总的方案数减去经过障碍的方案数。假设障碍的位置是 $(x, y)$，终点是 $(h, w)$，经过障碍的方案数就是 $C_{x+y-2}^{x-1} * C_{h-x+w-y}^{h-x}$（乘法原理）。

有多个障碍怎么办呢？联想到容斥原理，将总的方案数减去至少经过 $1$ 个障碍的，加上 $2$ 个的，减去 $3$ 个的……但这样做复杂度很高。

可以设 $dp_i$ 表示从 $(1, 1)$ 到 $(x_i, y_i)$ 且中途不经过其它障碍的方案数。令终点 $(h, w)$ 为第 $n+1$ 个障碍，求的答案就是 $dp_{n+1}$。

依然是正难则反，得出方程：$dp_i = C_{x_i+y_i-2}^{x_i-1}-\sum_{j=1}^{i-1}dp_j*C_{x_i-x_j+y_i-y_j}^{x_i-x_j}$。这样就可以巧妙地利用前面的状态不重不漏地计数了。

要写逆元。

 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
32
33
34
35
36
37
38
39
40

#include <bits/stdc++.h>
using namespace std;
using LL = long long;
const LL p = 1000000007;
int h, w, n;
LL dp[3005], fac[300010], inv[300010], invf[300010];
LL C(int n, int k) {
 if (n == k || k == 0) return 1;
 return fac[n] * invf[k] % p * invf[n - k] % p;
}
LL calc(int x1, int y1, int x2, int y2) { return C(x2 + y2 - x1 - y1, x2 - x1); }
struct za {
 int r, c;
 bool operator<(const za& b) {
 if (r != b.r) return r < b.r;
 return c < b.c;
 }
} a[3005];
int main() {
 ios::sync_with_stdio(false);
 inv[1] = invf[1] = fac[1] = 1;
 for (int i = 2; i <= 300000; i++) {
 inv[i] = (p - p / i) * inv[p % i] % p;
 invf[i] = invf[i - 1] * inv[i] % p;
 fac[i] = (fac[i - 1] * i) % p;
 }
 cin >> h >> w >> n;
 for (int i = 1; i <= n; i++) cin >> a[i].r >> a[i].c;
 sort(a + 1, a + 1 + n);
 a[++n] = {h, w};
 for (int i = 1; i <= n; i++) {
 dp[i] = calc(1, 1, a[i].r, a[i].c);
 for (int j = 1; j < i; j++) {
 if (a[j].r <= a[i].r && a[j].c <= a[i].c)
 dp[i] = (dp[i] - dp[j] * calc(a[j].r, a[j].c, a[i].r, a[i].c) % p + p) % p;
 }
 }
 cout << dp[n] << endl;
 return 0;
}

Orz：AT4546 题解 - GaryH

2023-05-05 08:00:00

我也不知道这是从哪本书上抠来的？

整除

定义 1：如果 $a$ 和 $b$ 为整数且 $a \ne 0$，我们说 $a$ 整除 $b$ 是指存在整数 $c$ 使得 $b=ac$。如果 $a$ 整除 $b$，我们还称 $a$ 是 $b$ 的一个因子，且称 $b$ 是 $a$ 的倍数。

如果 $a$ 整除 $b$，则将其记为 $a \mid b$，如果 $a$ 不能整除 $b$，则记其为 $a \nmid b$。

定理 1：如果 $a, b$ 和 $c$ 是整数，且 $a \mid b, b \mid c$，则 $a \mid c$。

定理 2：如果 $a, b, m$ 和 $n$ 为整数，且 $c \mid a, c \mid b$，则 $c \mid (ma+nb)$。

定理 3：带余除法 如果 $a$ 和 $b$ 是整数且 $b \gt 0$，则存在唯一的整数 $q$ 和 $r$，使得$a = bq + r, 0 ≤ r < b$。

最大公因子及其性质

定义 2：不全为零的整数 $a$ 和 $b$ 的最大公因子是指能够同时整除 $a$ 和 $b$ 的最大整数。

定义 3：设 $a,b$ 均为非零整数，如果 $a$ 和 $b$ 最大公因子 $(a,b)=1$，则称 $a$ 与 $b$ 互素。

定理 4：$a，b$ 是整数，且 $(a, b)=d$，那么 $(a/d, b/d)=1$。（换言之，$a/d$ 与 $b/d$ 互素）

推论 1：如果 $a, b$ 为整数，且 $b\ne 0$，则 $a/b=p/q$，其中 $p, q$ 为整数，且 $(p,q)=1, q≠0$。

定理 5：令 $a, b, c$ 是整数，那么 $(a+cb, b) = (a, b)$

证明：令 $a, b, c$ 是整数，证明 $a, b$ 的公因子与 $a+cb, b$ 的公因子相同，即证明 $(a+cb, b)=(a, b)$。

令 $e$ 是 $a, b$ 的公因子，由定理 2 可知 $e \mid (a+cb)$，所以 $e$ 是 $a+cb$ 和 $b$ 的公因子。

如果 $f$ 是 $a+cb$ 和 $b$ 的公因子，由定理 2 可知 $f$ 整除 $(a+cb)-cb=a$，所以 $f$ 是 $a, b$ 的公因子，因此 $(a+cb, b)=(a,b)$。

定义 4：线性组合 如果 $a,b$ 是整数，那么它们的线性组合具有形式 $ma+nb$，其中 $m,n$ 都是整数。

定理 6：两个不全为零的整数 $a, b$ 的最大公因子是 $a, b$ 的线性组合中最小的正整数。

证明：令 $d$ 是 $a,b$ 的线性组合中最小的正整数，$d = ma + nb$，其中 $m,n$ 是整数，我们将证明 $d\mid a, d\mid b$。

由带余除法，得到 $a=dq+r, 0\le r\lt d$。

由 $a=dq+r $和 $d=ma+nb$，得到 $r=a-dq=a-q(ma+nb)=(1-qm)a-qnb$。

这就证明了整数 $r$ 是 $a,b$ 的线性组合。因为 $0 \le r \lt d$，而 $d$ 是 $a,b$ 的线性组合中最小的正整数，

于是我们得到 $r=0$（如果 $r$ 不是等于 $0$，那意味着 $r$ 才是所有线性组合中最小的正整数，这与 $d$ 是所有线性组合中最小的正整数矛盾），因此 $d\mid a$，同理可得，$d\mid b$。

我们证明了 $a,b$ 的线性组合中最小的正整数 $d$ 是 $a,b$ 的公因子，剩下要证的是它是 $a,b$ 的最大公因子，为此只需证明 $a,b$ 所有的公因子都能整除 $d$。

由于 $d = ma + nb$，因此如果 $c \mid a$ 且 $c \mid b$，那么由定理 2 有 $c \mid d$，因此 $d \gt c$，这就完成了证明。

定义 5：令 $a_1,a_2,…,a_n$ 是不全为零的整数，这些整数的公因子中最大的整数就是最大公因子。$a_1,a_2,…,a_n$ 的最大公因子记为 $(a_1, a_2, …, a_n)$。

引理 1：如果 $a_1,a_2,…,a_n$ 是不全为零的整数，那么 $(a_1, a_2, …, a_{n-1}, a_n) = (a_1, a_2, …, (a_{n-1}, a_n))$。

欧几里得算法

辗转相除法求 gcd，即 $\gcd(a, b) = \gcd(b, a\bmod b)$。

证明 1：由定理 3 带余除法，存在整数 $q, r$ 使得 $a = bq + r, 0 \le r \lt b$, 得到 $r = a - bq$。由定理 5 得，$\gcd(a,b) = \gcd(a-bq,b) = \gcd(r,b) = \gcd(a%b,b) = \gcd(b,a%b)$。

证明 2：令 $d=(a,b)$，证明 $d \mid (b,a\bmod b)$，再反证 $(b,a\bmod b) \gt d$ 是不可能的。

时间复杂度的证明：

假设 $a\gt b$，分两种情况：

1. $b \lt a/2$， 经过一次辗转得到 $(b,a\bmod b)$，$\max(a,b)$ 至少缩小一半。
2. $b \ge a/2$，经过两次辗转得到彼时的 $\max(a,b)$ 至少缩小一半。

综上所述，最多 $2\log(n)$ 次辗转算法结束。

裴蜀定理

如果 $a$ 与 $b$ 均为整数，则存在整数 $x$ 和 $y$ 满足 $ax + by = (a,b)$。

由定理 6 容易证明正确性。

下面用扩展欧几里得算法（exgcd）求出满足 $ax + by = (a,b)$ 的一个特解。

例如：$a = 99, b = 78$， 令 $d =(a,b) = (99,78) = 3$，求 $99x + 78y = 3$ 的一个特解。

在调用 exgcd 的时候，从后往前依次构造出相应的解。

在用欧几里得算法求 $(99,78)$ 的时候，是要先求 $(78,21)$。

假如已经求出 $78x + 21y = 3$ 的一个解 ${x_0,y_0} = {3,-11}$，即 $78x_0 + 21y_0 = 3$。

那么可以由 $78x_0 + 21y_0 = 3$，构造出 $99x + 78y = 3$ 的一个特解。

因 $a=99, b=78, a\bmod b=21$, 因此 $78x_0 + 21y_0 = 3$，可以写成：$bx_0 + (a\bmod b)y_0 = 3$，即 $bx_0+(a-\frac{a}{b}b)y_0=3$，即 $ay_0+b(x_0-\frac{a}{b}y_0)=3$，即 $99y_0+78(x_0-\frac{99}{78}y_0)=3$。

那么只需要令 $x = y_0 = -11, y = x_0 - \frac{99}{78}y_0=14$，就可以得到 $99x + 78y = 3$ 的一个特解，即 ${-11, 14}$ 是 $99x+78y=3$ 的一个特解。

也就是说，在用欧几里得算法求 $(78,21)$ 的时候，若能返回 ${x_0,y_0}$ 使得满足 $78x_0 + 21y_0 = 3$，那么就能构造出一个特解 ${x,y}$ 满足 $99x + 78y = 3$ 的一个特解。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

void exgcd(int a, int b, int &d, int &x, int &y) {
 if (!b) {
 d = a;
 x = 1;
 y = 0;
 return;
 }
 exgcd(b, a % b, d, x, y);
 int t = x;
 x = y;
 y = t - (a / b) * y;
 return;
}

注意：若 $a\lt 0$ 且 $b\ge 0$ 那么在求 $ax+by=(a,b)$ 的时候，可以先求出 $|a|x+by=(|a|,b)$ 的一组解 ${x_0,y_0}$，然后 ${-x_0,y_0}$ 就是原方程的一个解。

若 $a\ge 0$ 且 $b\lt 0$ 的同理。若 $a \lt 0$ 且 $b \lt 0$ 的也同理。

定理 7：如果 $a,b$ 是正整数，那么所有 $a,b$ 的线性组合构成的集合与所有 $(a,b)$ 的倍数构成的集合相同。

证明：假设 $d = (a,b)$。

1. 首先证明每个 $a,b$ 的线性组合是 $d$ 的倍数。首先注意到由最大公因子的定义，有 $d\mid a$ 且 $d\mid b$，每个 $a,b$ 的线性组合具有形式 $ma+nb$，其中 $m,n$ 是整数。

由定理 2，只要 $m,n$ 是整数，$d$ 就整除 $ma+nb$，因此，$ma+nb$ 是 $d$ 的倍数。

2. 现在证明每一个 $d$ 的倍数也是 $(a，b)$ 的线性组合。

由定理 6，存在整数 $r,s$ 使得 $(a,b) = ra + sb$。而 $d$ 的倍数具有形式 $jd$，其中 $j$ 是整数。

在方程 $d = ra + sb$ 的两边同时乘以 $j$，我们得到 $jd = (jr)a + (js)b$。

因此，每个 $d$ 的倍数是 $(a,b)$ 的线性组合。

推论 2：整数 $a$ 与 $b$ 互素当且仅当存在整数 $m$ 和 $n$ 使得 $ma + nb = 1$。

证明：如果 $a，b$ 互素，那么 $(a,b)=1$，由定理 6 可知，$1$ 是 $a$ 和 $b$ 的线性组合的最小正整数，于是存在整数 $m,n$ 使得 $ma + nb = 1$。

反之，如果有整数 $m$ 和 $n$ 使得 $ma + nb = 1$，则由定理 6 可得 $(a,b)=1$，这是由于 $a,b$ 不同为 $0$ 且 $1$ 显然是 $a,b$ 的线性组合中的最小正整数。

二元一次不定方程

引理 2：如果 $a, b, c$ 是正整数，满足 $(a, b) = 1, a \mid bc$，则 $a \mid c$。

证明：由于 $(a, b)=1$，存在整数 $x$ 和 $y$ 使得 $ax+by=1$。等式两边同时乘以$c$，得 $acx+bcy=c$。

根据定理 2，$a$ 整除 $(cx)a + y(bc)$，这是因为这是 $a$ 和 $bc$ 的线性组合，而它们都可以被 $a$ 整除。因此，$a \mid c$。

定理 8：设 $a,b$ 是整数且 $d=(a,b)$。如果 $d \nmid c$，那么方程 $ax+by=c$ 没有整数解，如果 $d \mid c$，那么存在无穷多个整数解。

另外，如果 $x = x_0, y = y_0$ 是方程的一个特解，那么所有的解可以表示为：$x = x_0 + (b/d)n， y = y_0 - (a/d)n$，其中 $n$ 是整数。

证明：由定理 7 可知，$ax+by$ 的结果是 $d$ 的倍数，因此如果 $d \nmid c$，那么方程 $ax+by=c$ 没有整数解。

如果 $d\mid c$，存在整数 $s, t$ 使得 $as+bt=d$。

因为 $d\mid c$，存在整数 $e$ 使得 $de = c, c = de = (as+bt)e = a(se)+b(te)$

因此 $x_0 = se, y_0 = te$ 是方程 $ax + by = c$ 的一个特解。

为了证明方程存在无穷多个解，令 $x = x_0 + (b/d)n, y = y_0 - (a/d)n$，其中 $n$ 是整数。

1. 证明任何一对整数 $(x_0+(b/d)n, y_0 - (a/d)n)$ 它是方程的解。

因为 $a(x_0+(b/d)n) + b(y_0 - (a/d)n) = ax_0 + by_0 + a(b/d)n - b(a/d)n = ax_0 + by_0$，而 $ax_0+by_0$ 是方程 $ax+by=c$ 的解，所以 $(x_0+(b/d)n, y_0 -(a/d)n)$ 就是方程的解。

2. 证明方程的任何一个解都具有 $(x_0 + (b/d)n, y_0 - (a/d)n)$ 这种形式。

假设整数 $x,y$ 满足 $ax+by=c$，又因为 $ax_0+by_0=c$，两式相减得到：$a(x-x_0)+b(y-y_0)=0$。

即 $a(x-x_0) = b(y_0-y)$，等式两边同时除以 $d$ 得到 $(a/d)(x-x_0) = (b/d)(y_0-y)$，根据定理 4，$a/d$ 与 $b/d$ 互质，再根据引理 2，$(a/d) \mid (y_0-y)$，因此存在整数 $n$ 使得 $(a/d)n = y_0 - y$，于是得到 $y=y_0-(a/d)n$。

同理可得，$(b/d) \mid (x-x_0)$，因此存在整数 $n$ 使得 $(b/d)n = x - x_0$，于是得到 $x=x_0+(b/d)n$。

习题

1. Luogu-P5656 【模板】二元一次不定方程 (exgcd) Orz 离散小波变换°

下面设 $p = b/d, q = a/d$。

有正整数解：

• 解的个数：不断地将 $y_0$ 减去 $q$，$x_0$ 加上 $p$ 就可以找到所有可行解，个数为 $\lfloor (y-1)/q\rfloor + 1$
• $x$ 的最小正整数值：exgcd 得到的 $x_0$
• $y$ 的最小正整数值：不断地将 $y_0$ 减去 $q$（因为 $x_0$ 最小时，$y_0$ 一定最大），答案为 $(y-1)\bmod q+1$（特别注意 $0$ 的情况）
• $x$ 的最大正整数值：不断将 $x_0$ 加上 $p$，答案为 $x+\lfloor (y-1)/q\rfloor * p$
• $y$ 的最大正整数值：$x_0$ 为最小正整数时，$y_0$ 就是最大值

无正整数解：

• $x$ 的最小正整数值：exgcd 得到的 $x_0$
• $y$ 的最小正整数值：当前的 $y_0 \le 0$，需要执行构造 $x_0$ 的方法，即 $y_0 + q * \lceil (1.0-y)/q \rceil$

 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
32
33
34
35
36
37
38
39
40

#include <bits/stdc++.h>
using namespace std;
#define LL long long
void exgcd(LL a, LL b, LL &d, LL &x, LL &y) {
 if (b == 0) {
 d = a;
 x = 1;
 y = 0;
 return;
 }
 exgcd(b, a%b, d, x, y);
 LL t = x;
 x = y;
 y = t-(a/b)*y;
}
int main()
{
 int t;
 cin >> t;
 while (t--) {
 LL a, b, c, x, y, d;
 cin >> a >> b >> c;
 exgcd(a, b, d, x, y);
 if (c%d != 0) {
 puts("-1");
 } else {
 x *= c/d;
 y *= c/d;
 LL p = b/d, q = a/d, k;
 if (x < 0) { k = ceil((1.0-x)/p); x += p*k; y -= q*k; }
 if (x >= 0) { k = (x-1)/p; x -= p*k; y += q*k; }
 if (y > 0) {
 cout << (y-1)/q+1 << " " << x << " " << (y-1)%q+1 << " " << x+(y-1)/q*p << " " << y << endl;
 } else {
 cout << x << " " << y+q*(LL)ceil((1.0-y)/q) << endl;
 }
 }
 }
 return 0;
}

2. 多元一次不定方程的一组解：求 $a_1x_1 + a_2x_2 + a_3x_3 + … + a_nx_n = c$ 的一组整数解，如无整数解输出 $-1$。

先考虑三元一次不定方程 $a_1x_1 + a_2x_2 + a_3+x_3 = c$。可以用 exgcd 解二元一次方程 $a_1x_1 + a_2x_2 = (a_1, a_2)$，设 $d = (a_1, a_2)$，由定理 8 知道 $d$ 乘任何整数时该方程都有整数解，就可以把 $d$ 当做原方程的一个系数，转而求解 $dt+a_3x_3 = c$。上一步的方程变为 $t(a_1x_1 + a_2x_2) = td$，于是将 $t$ 乘上前面的 $x_1, x_2$ 就得到最终答案。

多元同理。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

int main()
{
 exgcd(a[1], a[2], ta[2], ans[1], ans[2]);
 int tm = 1;
 for (int i = 3; i <= n; i++) {
 LL tmp;
 exgcd(ta[i-1], a[i], ta[i], tmp, ans[i]);
 if (i == n && c%ta[i]) { cout << -1 << endl; return 0; }
 for (int j = 1; j < i; j++) ans[j] *= tmp;
 }
 tm = c/ta[n];
 for (int i = 1; i <= n; i++) {
 cout << ans[i]*tm <<" ";
 }
 return 0;
}

3. Luogu-P3986 斐波那契数列：求 $f(i)x+f(i+1)y=k$ 的正整数解个数（$f$ 表示斐波那契数列）

显然 $f(i)+f(i+1)\gt k$ 时就不用继续下去了，因此方程的个数是有限的，可以枚举。

然后题目就变成了求 $f_ia+f_{i+1}b=k$ 的正整数解的个数。

另外，有没有可能同样的 $(a, b)$ 出现了两次呢？不可能。否则就需要满足 $af_i+bf_{i+1}=af_j+bf_{j+1}, i \ne j$，然而斐波那契数列任两项不相等，以上式子不成立。

 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

int main()
{
 LL n;
 cin >> n;
 fib[0] = fib[1] = 1;
 for (int i = 2; i <= 50; i++) fib[i] = fib[i-1]+fib[i-2];
 LL ans = 0;
 for (int i = 1; fib[i] < n; i++) {
 LL d, x, y, k;
 exgcd(fib[i-1], fib[i], d, x, y);
 LL p = fib[i]/d, q = fib[i-1]/d;
 x = x*(n/d), y = y*(n/d);
 if (x <= 0) {
 k = ceil((1.0-x)/p);
 x += k*p;
 y -= k*q;
 }
 if (x >= 0) {
 k = (x-1)/p;
 x -= k*p;
 y += k*q;
 }
 if (x <= 0 || y <= 0) continue;
 ans = ((y-1)/q+1+ans)%MOD;
 }
 cout << ans << endl;
 return 0;
}

同余

同余概述

定义 6：设 $m$ 是正整数，若 $a$ 和 $b$ 是整数，且 $m\mid (a-b)$，则称 $a$ 和 $b$ 模 $m$ 同余。

若 $a$ 和 $b$ 模 $m$ 同余，则记 $a\equiv b(mod\ m)$。

若 $m \nmid (a-b)$，则记 $a\not\equiv b (mod\ m)$，并称 $a$ 模 $m$ 不同余于 $b$。

整数 $m$ 称为同余的模。

定理 9：若 $a$ 和 $b$ 是整数，则 $a\equiv b(mod\ m)$ 当且仅当存在整数 $k$，使得 $a=b+km$。

证明：若 $a\equiv b(mod\ m)$，则 $m\mid (a-b)$，这说明存在整数 $k$，使得 $km=a-b$，所以 $a=b+km$。

反过来，若存在整数 $k$ 使得 $a=b+km$，则 $km=a-b$。于是，$m\mid (a-b)$，因而 $a\equiv b(mod\ m)$。

例：$19\equiv -2(mod\ 7)$ 和 $19=-2+3·7$。

定理 10：设 $m$ 是正整数，模 $m$ 的同余满足下面的性质:

1. 自反性。若 $a$ 是整数，则 $a\equiv a(mod\ m)$。
2. 对称性。若 $a$ 和 $b$ 是整数，且 $a\equiv b(mod\ m)$，则 $b\equiv a(mod\ m)$。
3. 传递性。若 $a,b,c$ 是整数，且 $a\equiv b(mod\ m)$ 和 $b\equiv c(mod\ m)$，则 $a\equiv c(mod\ m)$。

证明：

1. 因为 $m\mid(a-a)=0$，所以 $a=a(mod\ m)$。
2. 若 $a\equiv b(mod\ m)$，则 $m\mid(a-b)$。从而存在整数 $k$，使得 $km=a-b$。这说明 $(-k)m=b-a$，即 $m\mid (b-a)$。因此，$b\equiv a(mod\ m)$。
3. 若 $a\equiv b(mod\ m)$，且 $b\equiv c(mod\ m)$，则有 $m\mid (a-b)$ 和 $m\mid (b-c)$。从而存在整数 $k$ 和 $l$，使得 $km=a-b,lm=b-c$。于是，$a-c=(a-b)+(b-c)=km+lm=(k+l)m$。因此，$m\mid(a-c), a\equiv c(mod\ m)$。

由定理 10 可见，整数的集合被分成 $m$ 个不同的集合，这些集合称为模 $m$ 剩余类（同余类），每个同余类中的任意两个整数都是模 $m$ 同余的。

注意，当 $m=2$ 时，正好整数分成奇、偶两类.

如果你对集合上的关系比较熟悉，那么定理 10 表明对正整数 $m$ 的模 $m$ 同余是一种等价关系，并且每一个模 $m$ 同余类即是由此种等价关系所定义的等价类。

例模 $4$ 的四个同余类是：

$$ …\equiv 8\equiv -4\equiv 0\equiv 4\equiv 8\equiv …(mod\ 4) \\ …\equiv -7\equiv -3\equiv 1\equiv 5\equiv 9\equiv …(mod\ 4) \\ …\equiv -6\equiv -2\equiv 2\equiv 6\equiv 10\equiv …(mod\ 4) \\ …\equiv -5\equiv -1\equiv 3\equiv 7\equiv 11\equiv …(mod\ 4) \\ $$

设 $m$ 是正整数，给定整数 $a$，由带余除法有 $a = bm + r$，其中 $0\le r \lt m-1$，称 $r$ 为 $a$ 的模 $m$ 最小非负剩余，是 $a$ 模 $m$ 的结果，类似地，当 $m$ 不整除 $a$ 时，称 $r$ 为 $a$ 的模 $m$ 最小正剩余。

$mod\ m$ 实际上是从整数集到集合 ${0,1,2,…,m-1}$ 的函数。

定理 11：若 $a$ 与 $b$ 为整数，$m$ 为正整数，则 $a\equiv b(mod\ m)$ 当且仅当 $a \bmod m = b \bmod m$。

证明：

1. 若 $a\equiv b(mod\ m)$ 则 $a\bmod m = b \bmod m$。

由带余除法 $a=pm+ra , b=qm+rb$，因 $a\equiv b(mod\ m)$，则 $m\mid a-b$，则 $m\mid (p-q)m+(ra-rb)$，则存在整数 $x$，满足 $xm = (p-q)m+(ra-rb)$，则 $ra-rb = (x-p+q)m$，则 $ra-rb\mid m$。

因 $0\le ra\lt m, 0\le rb\lt m$, 故 $-(m-1) \le ra-rb \lt m$，故 $ra-rb$ 只能是 $0$ 才能满足 $ra-rb\mid m$。

则 $ra = rb$，则 $ra=a \bmod m, rb=b \bmod m$，则 $a \bmod m = b \bmod m$。

2. 若 $a \bmod m = b \bmod m$ 则 $a\equiv b(mod\ m)$。

由带余除法 $a=pm+ra , b=qm+rb$，若 $a \bmod m = b \bmod m$，则 $ra=rb$。

因此 $a-b = (pm+ra)-(qm+rb)=(p-q)m + (ra-rb) = (p-q)m$。因此 $m\mid a-b$，故 $a\equiv b(mod\ m)$。

因此，每个整数都和 $0,1,…,m-1$（也就是 $a$ 被 $m$ 除所得的余数）中的一个模 $m$ 同余。因为 $0,1,…,m-1$ 中的任何两个都不是模 $m$ 同余的，所以有 $m$ 个整数使得每个整数都恰与这 $m$ 个整数中的一个同余。

定理 12：若 $a,b,c,m$ 是整数，$m\gt 0$，且 $a\equiv b(mod\ m)$，则

1. $a+c\equiv b+c(mod\ m)$
2. $a-c\equiv b-c(mod\ m)$
3. $ac\equiv bc(mod\ m)$

定理 13：若 $a,b,c,m$ 是整数，$m\gt 0, d=(c,m)$，且有 $ac\equiv bc(mod\ m)$，则 $a\equiv b(mod\ m/d)$。

证明：若 $ac\equiv bc(mod\ m)$，所以存在整数 $k$ 满足 $c(a-b)=km$，两边同时除以 $d$，得到：$(c/d)(a-b)=k(m/d)$，因为 $(m/d,c/d)=1$，根据引理 2，$m/d\mid a-b$，因此 $a\equiv b(mod\ m/d)$。

下面的推论是定理 13 的特殊情形，经常用到，它使得我们能够在模 $m$ 同余式中消去与模 $m$ 互素的数。

推论 3：若 $a,b,c,m$ 是整数，$m\gt 0,(c，m)=1$，且有 $ac\equiv bc(mod\ m)$，则 $a\equiv b(mod\ m)$。

定理 14 ：若 $a,b,c,m$ 是整数，$m\gt 0,a\equiv b(mod\ m)$，且 $c\equiv d(mod\ m)$，则:

1. $a+c\equiv b+d(mod\ m)$
2. $a-c\equiv b-d(mod\ m)$
3. $ac\equiv bd(mod\ m)$

证明：

因为 $a\equiv b(mod\ m)$ 且 $c\equiv d(mod\ m)$，我们有 $m\mid (a-b)$ 与 $m\mid (c-d)$。因此，存在整数 $k$ 与 $l$ 使得 $km=a-b,lm=c-d$。

为证 1，注意到 $(a+c)-(b+d)=(a-b)+(c-d)=km+lm=(k+l)m$.因此 $m\mid (a+c)-(b+d)$，即 $a+c\equiv b+d(mod\ m)$。

为证 2，注意到 $(a-c)-(b-d)=(a-b)-(c-d)=km-lm=(k-1)m$，因此 $m\mid (a-c)-(b-d)$，即 $a-c\equiv b-d(mod\ m)$.

为证 3，注意到 $ac-bd=ac-bd+bc-bd=c(a-b)+b(c-d)=ckm+blm=m(ck+ bl)$。因此，$m\mid ac-bd$，即 $ac\equiv bd(mod\ m)$。

定义 7：一个模 $m$ 完全剩余系是一个整数的集合，使得每个整数恰和此集合中的一个元素模 $m$ 同余。

例如：整数 $0,1,…,m-1$ 的集合是模 $m$ 完全剩余系，称为模 $m$ 最小非负剩余的集合。

下面的引理帮助我们判定一个 $m$ 元集合是否为模 $m$ 的完全剩余系.

引理 3：$m$ 个模 $m$ 不同余的整数的集合构成一个模 $m$ 的完全剩余系。

证明：

假设 $m$ 个模 $m$ 不同余的整数集合不是模 $m$ 完全剩余系，这说明，至少有一个整数 $a$ 不同余于此集合中的任一整数。

所以此集合中的整数都模 $m$ 不同余于 $a$ 被 $m$ 除所得的余数，从而，整数被 $m$ 除得的不同剩余至多有 $m-1$ 个。

由鸽笼原理，此集合中至少有两个整数有相同的模 $m$ 剩余。

这不可能，因为这些整数均模 $m$ 不同余，因此，$m$ 个模 $m$ 不同余的整数的集合构成一个模 $m$ 的完全剩余系。

定理 15：若 $r_1,r_2,…,r_m$ 是一个模 $m$ 的完全剩余系，且正整数 $a$ 满足 $(a,m)=1$，则对任何整数 $b$，$ar_1+b, ar_2+b,…,ar_m+b$ 都是模 $m$ 的完全剩余系。

证明：首先来证整数 $ar_1+b, ar_2+b,…,ar_m+b$ 中的任何两个都模 $m$ 不同余。

反证，若存在 $1\le j,k\le m$ 且 $j\ne k$ 且 $ar_j+b \equiv ar_k+b(mod\ m)$，则由定理 12.2 可知：$ar_j \equiv ar_k(mod\ m)$。因为 $(a，m)=1$，推论 3 表明 $r_j\equiv r_k(mod\ m)$，这与 $r_1,r_2,…,r_m$ 是一个模 $m$ 的完全剩余系相矛盾。

故 $ar_1+b, ar_2+b,…,ar_m+b$ 是 $m$ 个模 $m$​ 不同余的整数，由引理 3，命题得证。

定理 16：若 $a,b,k,m$ 是整数，$k\gt 0,m\gt 0$，且 $a\equiv b(mod\ m)$，则 $a^k\equiv b^k(mod\ m)$。

证明：因为 $a\equiv b(mod\ m)$，所以 $m\mid a-b$。

因为 $a^k-b^k =(a-b)(a^{k-1}+a^{k-2}b+…ab^{k-2}+b^{k-1})$ （可以参考资料：详聊如何理解a^n-b^n因式分解）

所以 $a-b \mid a^k-b^k$，根据 定理 1，$m \mid a^k-b^k$，即 $a^k\equiv b^k(mod\ m)$。

线性同余方程

设 $x$ 是未知整数，形如 $ax\equiv b(mod\ m)$ 的同余式称为一元线性同余方程。

首先注意到，若 $x=x_0$ 是同余方程 $ax\equiv b(mod\ m)$ 的一个解，且 $x_1\equiv x_0(mod\ m)$，则 $ax_1\equiv ax_0 \equiv b(mod\ m)$，所以 $x_1$ 也是一个解。

因此，若一个模 $m$ 同余类的某个元素是解，则此同余类的所有元素都是解。

于是，我们会问模 $m$ 的 $m$ 个同余类中有多少个是给出方程的解，这相当于问方程有多少个模 $m$ 不同余的解。

定理 17：设 $a,b,m$ 是整数，$m\gt 0,(a,m)=d$。

若 $d\nmid b$，则 $ax\equiv b(mod\ m)$ 无解。

若 $d\mid b$，则 $ax\equiv b(mod\ m)$ 恰有 $d$ 个模 $m$ 不同余的解。

证明：根据定理 9，线性同余方程 $ax\equiv b(mod\ m)$ 可以写成二元线性不定方程 $ax+my=b$。其中 $x, y$ 是未知数。整数 $x$ 是 $ax\equiv b(mod\ m)$ 的解当且仅当存在整数 $y$ 使得 $ax+my=b$。

由定理 8 可知，若 $d\nmid b$，则无解，而 $d\mid b$ 时，$ax-my=b$ 有无穷多解：$x = x_0 + (m/d)t, y = y - (a/d)t$，

其中 $x=x_0$ 和 $y=y_0$ 是方程的特解，上述 $x$ 的值 $x=x_0+(m/d)t$ 是线性同余方程的解，有无穷多这样的解。

为确定有多少不同余的解，我们来找两个解 $x_1=x_0+(m/d)t1$ 和 $x_2=x_0+(m/d)t2$ 模 $m$ 同余的条件，若这两个解同余，则 $(x0+(m/d)t1)\equiv (x0+(m/d)t2) (mod\ m)$。

根据 定理 12，两边减去 $x_0$，有 $(m/d)t1\equiv (m/d)t2(mod\ m)$。

因为 $(m/d)\mid m$，所以 $(m,m/d)=m/d$，再由 定理 13 得，$t1\equiv t2 (mod\ d)$。

这表明不同余的解的一个完全集合可以通过取 $x=x_0+(m/d)t$ 得到，其中 $t$ 取遍模 $d$ 的完全剩余系，一个这样的集合可由 $x=x_0+(m/d)t$ 给出，其中 $t=0,1,2,…,d-1$。

推论 4：若 $a$ 和 $m\gt 0$ 互素，且 $b$ 是整数，则线性同余方程 $ax\equiv b(mod\ m)$ 有模 $m$ 的唯一解。

证明：

因为 $(a,m)=1$，所以 $(a,m)\mid b$。因此，由 定理 17，线性同余方程 $ax\equiv b(mod\ m)$ 恰有 $(a,m)=1$ 个模 $m$ 不同余的解。

例：为求出 $9x\equiv 12(mod\ 15)$ 的所有解，首先注意到因为 $(9,15)=3$ 且 $3\mid 12$，所以恰有三个不同余的解，我们可以通过先找到一个特解，再加上 $15/3=5$ 的适当倍数来求得所有的解。

为求特解，我们考虑二元线性不定方程 $9x-15y=12$。由扩展欧几里得算法得到：$9x-15y=12$ 的一个特解是 $x_0=8$ 和 $y_0=4$。

由定理 17 的证明可知，三个不同余的解由 $x= x_0\equiv 8(mod\ 15),x= x_0+5\equiv 13(mod\ 15)$ 和 $x= x_0+5*2\equiv 18\equiv 3(mod\ 15)$ 给出。

习题

Luogu-P1516 青蛙的约会 Orz 题解 P1516 【青蛙的约会】 - FlashHu

如果两蛙相遇，那么他们的初始坐标差 $x-y$ 和跳的距离 $(n-m)t$ 之差应该模纬度线总长 $l$ 同余，$(n-m)t\equiv x-y(mod\ l)$。转化成不定方程的形式：$(n-m)t+kl=x-y$，并求最小正整数解。设 $a=n-m,b=l,c=x-y$，可以写成 $ax+by = c$，通过 exgcd 可以求出 x 的一个特解。

细节问题，因为 gcd 只对非负整数有意义，如果 $a\lt 0$ 时等式两边要同时取负，$a,c$ 变成相反数（相当于把两个蛙交换了一下），$b$ 是正数所以不能变。

这里求最小正整数解时用了模的方法来处理，值得细品 (x0%p+p)%p。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

int main()
{
 cin >> x >> y >> m >> n >> l;
 LL a = n-m, b = l, c = x-y, d, x0, y0;
 if (a < 0) {
 a = -a;
 c = -c;
 }
 exgcd(a, b, d, x0, y0);
 if (c % d == 0) {
 x0 *= c/d;
 y0 *= c/d;
 LL p = b/d;
 cout << (x0%p+p)%p << endl;
 } else {
 cout << "Impossible\n";
 }
 return 0;
}

POJ-2115 – C Looooops

可将题意转化为 $A+Ct\equiv B(mod\ 2^K)$，把它转化成方程：$A+Ct-B=2^Kz$，即 $Ct+2^Kz=B-A$，用 exgcd 解这个方程并求出 $t$ 的最小正整数解即为答案。

注意 1LL<<K，不开 long long 差点怀疑自己做错了。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

while (cin >> A >> B >> C >> K && !(A==0&&B==0&&C==0&&K==0)) {
 K = 1ll<<K;
 LL a = C, b = K, d, x, y;
 exgcd(a, b, d, x, y);
 if (((B-A+K)%K)%d) {
 cout <<"FOREVER\n";
 } else {
 x *= ((B-A+K)%K)/d;
 LL p = b/d;
 x = (x%p+p)%p;
 cout << x << endl;
 }
}

模的逆

现在考虑特殊形式的同余方程 $ax\equiv 1(mod\ m)$。

由 定理 17，此方程有解当且仅当 $(a,m)=1$，于是其所有的解都模 $m$ 同余。

定义 8：给定整数 $a$，且满足 $(a,m)=1$，称 $ax\equiv 1(mod\ m)$ 的一个解为 $a$ 模 $m$ 的逆。

例：因为 $7x\equiv 1(mod\ 31)$ 的解满足 $x\equiv 9(mod\ 31)$，所以 $9$ 和所有与 $9$ 模 $31$ 同余的整数都是 $7$ 模 $31$ 的逆。类似地，因为 $9·7\equiv 1(mod\ 31)$，所以 $7$ 是 $9$ 模 $31$ 的逆。

当我们有 $a$ 模 $m$ 的一个逆时，可以用它来解形如 $ax\equiv b(mod\ m)$ 的任何同余方程。

为看清这一点，令 $\bar{a}$ 是 $a$ 的模 $m$ 的一个逆，所以 $a\bar{a}\equiv 1(mod\ m)$。

于是，若 $ax\equiv b(mod\ m)$，则根据 定理 12 将同余方程两边同时乘以 $\bar{a}$，得到 $\bar{a}(ax)\equiv \bar{a}b(mod\ m)$，所以 $x=\bar{a}b(mod\ m)$。

例 为求出 $7x\equiv 22(mod\ 31)$ 的所有解，我们在此方程两边同时乘以 $9$（这是 $7$ 模 $31$ 的一个逆），得 $9·7x\equiv 9·22(mod\ 31)$。因此，$x\equiv 198=12(mod\ 31)$。

求逆的三种算法：

1. 扩展欧几里得算法解同余方程，求单个逆。$ax\equiv 1(mod\ m)$ 等价于求二元线性不定方程的解：$ax+my=1$，其中 $(a,m)=1$ 是有逆的前提。只需要用扩展欧几里得算法求即可。

2. 费马小定理，求单个逆

定理 18（费马小定理）：设 $p$ 是一个素数，$a$ 是一个正整数且 $p\nmid a$，则 $a^{p-1}\equiv 1(mod\ p)$。

证明：考虑 $p-1$ 个整数 $a,2a,…,(p-1)a$，它们都不能被 $p$ 整除。因为若 $p \mid ja$, 那么因 $p \nmid a$，则由 引理 2 知 $p \mid j$，但 $1\le j\le p-1$，故这是不可能的。

现证明 $a,2a,…,(p-1)a$ 中任何两个整数模 $p$ 不同余。

为了证明这一点，设 $ja\equiv ka(mod\ p)$，其中 $1\le j\lt k\le p-1$。

那么根据 推论 3，因 $(a,p)=1$，故 $j\equiv k(mod\ p)$，但这也是不可能的，因为 $j$ 和 $k$ 都是小于 $p-1$ 的正整数.

因为整数 $a,2a,…,(p-1)a$ 是 $p-1$ 个满足模 $p$ 均不同余于 $0$ 且任何两个都互不同余的整数组成的集合中的元素，故由 引理 3 可知 $a,2a,…,(p-1)a$ 模 $p$ 的最小正剩按照一定的顺序必定是整数 $1,2,…,p-1$。

由同余性，整数 $a,2a,…,(p-1)a$ 的乘积模 $p$ 同余于前 $p-1$ 个正整数的乘积，即：

$a·2a·3a···(p-1)a \equiv 1·2·3···(p-1) (mod\ p)$

因此 $a^{p-1}(p-1)! \equiv (p-1)! (mod\ p)$， 因 $(p,(p-1)!)=1$，根据 推论 3，消去 $(p-1)!$，得到 $a^{p-1}\equiv 1(mod\ p)$，证毕。

利用费马小定理，$a^{p-1}\equiv 1(mod\ p)$，则 $a·a^{p-2} \equiv 1(mod\ p)$，即 $a^{p-2}$ 是 $a$ 模 $p$ 的一个逆。

注意前提：$p$ 是一个素数，$a$ 是一个正整数且 $p \nmid a$。

通常可以用快速幂求 $a^{p-2}$。

3. 若 $p$ 是素数，$n\lt p$，线性递推求 $1$ 至 $n$ 在模 $p$ 意义下的乘法逆元

首先，$i=1$ 的逆是 $1$。下面求 $i\gt 1$ 的逆，用递推法。

1. 用带余除法 $p = k·i + r$，其中 $0\le r \lt i$，故 $k·i + r \equiv 0 (mod\ p)$
2. 在等式两边乘 $i^{-1}r^{-1}$，即 $(k·i + r)·i^{-1}r^{-1} \equiv 0 (mod\ p)$，即 $kr^{-1} + i^{-1} \equiv 0 (mod\ p)$
3. 移项得 $i^{-1}\equiv -kr^{-1} (mod\ p)$，即 $i^{-1}\equiv (-p/i)r^{-1} (mod\ p)$。若要避免负数，因 $pr^{-1} \equiv 0 (mod\ p)$，由定理 12 得，$pr^{-1} + (-p/i)r^{-1} \equiv (-p/i)r^{-1} (mod\ p)$ ，即 $(p-p/i)r^{-1} \equiv (-p/i)r^{-1} (mod\ p) $。

则 $i^{-1}\equiv (p-p/i)r^{-1} (mod\ p)$。

代码：

1
2
3
4

void inverse(LL n, LL p) {
 inv[1] = 1;
 for (int i = 2; i <= n; i++) inv[i] = (p-p/i)*inv[p%i]%p;
}

逆与除法取模

逆的一个重要应用是求除法的模。求 $(a/b) \bmod m$，即 $a$ 除以 $b$，然后对 $m$ 取模。

这里 $a$ 和 $b$ 都是很大的数，如 $a=n!$，容易溢出，导致取模出错。

用逆可以避免除法计算，设 $b$ 的逆元是 $b^{-1}$，有 $(a/b) \bmod m = ((a/b) \bmod m) ((bb^{-1}) \bmod m) = (a/b×bb^{-1}) \bmod m = (ab^{-1}) \bmod m$

经过上述推导，除法的模运算转换为乘法模运算，即

$(a/b) \bmod m = (ab^{-1}) \bmod m = (a \bmod m) (b^{-1} \bmod m) \bmod m$

习题

HDU-1576 A/B 因为 $\gcd(B,9973)=1$，可以用 exgcd 求逆元。

1
2
3
4
5
6
7
8
9

while (t--) {
 LL n, B;
 cin >> n >> B;
 LL a, b, d, x, y;
 exgcd(B, 9973, d, x, y);
 LL p = 9973;
 x = (x%p+p)%p;
 cout << n*x%9973 << endl;
}

更多资料

• 数学：Exgcd 与乘法逆元的漫步 - Aw顿顿の小窝 · Linear Expectation 线性期望 【数论】Exgcd/乘法逆元 - 题单

2023-04-30 08:00:00

Luogu-P4755 Beautiful Pair

题意

小 D 有个数列 ${a}$，当一个数对 $(i,j)$（$i \le j$）满足 $a_i$ 和 $a_j$ 的积不大于 $a_i, a_{i+1}, \ldots, a_j$ 中的最大值时，小 D 认为这个数对是美丽的。请你求出美丽的数对的数量。

$1\le n\le{10}^5$，$1\le a_i\le{10}^9$。

编程时的问题

• 对 ST 表不熟悉！
• 更 zz 的是，对 lower_bound 和 upper_bound 理解有问题，来复习一下小学知识：lower_bound 是找到“大于等于”的位置，upper_bound 是“大于”。写这道题的时候找小于某数的位置莫名其妙地用了 lower_bound，更没有 -1，完全是随手写的，半天也没察觉到这里有问题。

综上，我是 zz。

思路

考虑分治（据说这是套路），我们找出一个区间 $[l, r]$ 内的最大值位置 $mid$，然后统计所有跨过 $mid$ 的答案，再递归处理 $[l, mid-1], [mid+1, r]$。假设 $mid$ 左边的数是 $a_i$，右边的数是 $a_j$，根据题目得 $a_i * a_j \le a_{mid}$，即 $a_j \le \lfloor\frac{a_{mid}}{a_i}\rfloor$。那么我们枚举 $a_i$，然后用主席树统计右区间内小于 $\lfloor\frac{a_{mid}}{a_i}\rfloor$ 的数的个数。

注意每次要枚举左右区间中长度较小的那个，这样可以做到 $O(n\log^2n)$。否则会被卡成 $O(n^2\log 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
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84

#include <bits/stdc++.h>
using namespace std;
const int N = 100005, LOGN = 18;
int n, a[N], b[N], bn, f[N][LOGN], logn[N], root[N], cnt;
struct node {
 int sum, ls, rs;
} t[N*25];
void update(int &rt, int pre, int pos, int l, int r) {
 rt = ++cnt;
 t[rt] = t[pre];
 t[rt].sum++;
 if (l == r) return;
 int mid = (l+r)>>1;
 if (pos <= mid) update(t[rt].ls, t[pre].ls, pos, l, mid);
 else update(t[rt].rs, t[pre].rs, pos, mid+1, r);
}
int query(int rt, int pre, int x, int y, int l, int r) {
 if (l >= x && r <= y) return t[rt].sum-t[pre].sum;
 int mid = (l+r)>>1;
 int ret = 0;
 if (x <= mid) ret += query(t[rt].ls, t[pre].ls, x, y, l, mid);
 if (y > mid) ret += query(t[rt].rs, t[pre].rs, x, y, mid+1, r);
 return ret;
}
void buildst() {
 logn[1] = 0;
 logn[2] = 1;
 for (int i = 3; i < N; i++) logn[i] = logn[i/2]+1;
 for (int i = 1; i <= n; i++) f[i][0] = i;
 for (int t = 1; t < LOGN; t++) {
 for (int i = 1; i+(1<<t)-1 <= n; i++) {
 if (a[f[i][t-1]] > a[f[i+(1<<(t-1))][t-1]])
 f[i][t] = f[i][t-1];
 else
 f[i][t] = f[i+(1<<(t-1))][t-1];
 }
 }
}
int getmax(int l, int r) {
 int t = logn[r-l+1];
 // NOT f[l+(1<<(t-1))-1][t]
 if (a[f[l][t]] > a[f[r-(1<<t)+1][t]]) return f[l][t];
 return f[r-(1<<t)+1][t];
}
long long ans = 0;
void solve(int l, int r) {
 if (l > r) return;
 if (l == r) { ans += (a[l]==1); return; }
 int mid = getmax(l, r);
 if (mid-l+1 <= r-mid+1) { // 枚举左半区间
 for (int i = l; i <= mid; i++) {
 // 不能用lowerbound
 int idx = upper_bound(b+1, b+1+bn, a[mid]/a[i])-b-1;
 if (idx != 0)
 ans += query(root[r], root[mid-1], 1, idx, 1, bn);
 }
 } else { // 枚举右半区间
 for (int i = mid; i <= r; i++) {
 int idx = upper_bound(b+1, b+1+bn, a[mid]/a[i])-b-1;
 if (idx != 0)
 ans += query(root[mid], root[l-1], 1, idx, 1, bn);
 }
 }
 solve(l, mid-1);
 solve(mid+1, r);
}
int main()
{
 ios::sync_with_stdio(false);
 cin >> n;
 for (int i = 1; i <= n; i++)
 cin >> a[i], b[i] = a[i];
 sort(b+1, b+1+n);
 bn = unique(b+1, b+1+n)-b-1;
 root[0] = ++cnt;
 buildst();
 for (int i = 1; i <= n; i++) {
 int id = lower_bound(b+1, b+1+bn, a[i])-b;
 update(root[i], root[i-1], id, 1, bn);
 }
 solve(1, n);
 cout << ans << endl;
 return 0;
}

2023-04-28 08:00:00

消元法及高斯消元法思想

消元法是将方程组中的一方程的未知数用含有另一未知数的代数式表示，并将其带入到另一方程中，这就消去了一未知数，得到一解；或将方程组中的一方程倍乘某个常数加到另外一方程中去，也可达到消去一未知数的目的。消元法主要用于二元一次方程组的求解。

消元法理论的核心

消元法理论的核心主要如下：

• 两方程互换，解不变；
• 一方程乘以非零数 $k$，解不变；
• 一方程乘以数 $k$ 加上另一方程，解不变。

过程

解方程组：

$$ \begin{cases} 2x_1+x_2-x_3=8 \ -3x_1-x_2+2x_3=-11 \ -2x_1+x_2+2x_3=-3 \end{cases} $$

写成矩阵的形式为：

$$ \left[\begin{matrix} 2 & 1 & -1 \ -3 & -1 & 2 \ -2 & 1 & 2 \end{matrix} \middle| \begin{matrix} 8 \ -11 \ -3 \end{matrix} \right] $$

这种矩阵称为增广矩阵。所谓增广矩阵，即为方程组系数矩阵 $A$ 与常数列 $b$ 的并生成的新矩阵，即 $(A | b)$，增广矩阵行初等变换化为行最简形，即是利用了高斯消元法的思想理念，省略了变量而用变量的系数位置表示变量，增广矩阵中用竖线隔开了系数矩阵和常数列，代表了等于符号。

我们从上到下依次处理每一行，处理完第 $i$ 行后，让 $A_{ii}$ 非 $0$，而 $A_{ji}(j\gt i)$ 均为 $0$。过程如下。

$$ \left[\begin{matrix} 2 & 1 & -1 \ -3 & -1 & 2 \ -2 & 1 & 2 \end{matrix} \middle| \begin{matrix} 8 \ -11 \ -3 \end{matrix} \right] \

\Rightarrow

\left[\begin{matrix} -3 & -1 & 2 \ & 1/3 & 1/3 \ & 5/3 & 2/3 \end{matrix} \middle| \begin{matrix} -11 \ 2/3 \ 13/3 \end{matrix} \right] \

\Rightarrow

\left[\begin{matrix} -3 & -1 & 2 \ & 5/3 & 2/3 \ & & 1/5 \end{matrix} \middle| \begin{matrix} -11 \ 13/3 \ -1/5 \end{matrix} \right] $$

其中最后一个增广矩阵的系数部分是上三角阵（$0$ 用空白表示），而且主对角元 $a_{ii}$ 均非 $0$。该矩阵对应的方程组为：

$$ \begin{cases} -3x_1-x_2+2x_3=-11 \ 5x_2/3+2x_3/3=13/3 \ x_3/5=-1/5 \end{cases} $$

可得 $x_3=-1$，再代入前面两个方程可求解。

在消元部分中，假设正在处理第 $i$ 行，则首先需要找一个 $r\gt i$ 且绝对值最大的 $a_{ri}$，然后交换第 $r$ 行和第 $i$ 行。交换两个方程的位置不会对解产生影响，但可以减小计算误差。当 $A$ 可逆时，可以保证交换后 $a_{ii}$ 一定不等于 $0$。这种方法称为列主元法。

接下来进行所谓的“加减消元”。比如在上面的例子中的第一步，首先交换第一个和第二个方程的位置，然后用第一个方程 $(-3, 1, 2-11)$ 消去第二个方程 $(2, 1,-1, 8)$ 的第一列，方法是把第二个方程中的每个数都减去第一行对应元素的 $-2/3$ 倍。

一般情况下，如果要用第 $i$ 个方程来消去第 $k$ 个方程的第 $i$ 列，那么第 $k$ 行的所有元素 $A[k][j]$ 都应该减去 $A[i][j]$ 的 $A[k][i]/A[i][i]$ 倍。

下一个过程是回代。现在 $A$ 已经是一个上三角矩阵了，即第 $1,2,3…$ 行的最左边非 $0$ 元素分别在第 $1,2,3…$ 列。这样，最后一行实际上已经告诉我们 $x_n$ 的值了；接下来像前面说的那样不停地回代计算，最终会得到每个变量的唯一解。代码如下。

Luogu-P3389 【模板】高斯消元法

 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
32
33
34
35
36
37
38
39
40
41
42

#include <bits/stdc++.h>
using namespace std;
const int N = 106;
const double EPS = 1e-7;
int n;
double a[N][N];
bool gauss() {
 int r;
 for (int i = 0; i < n; i++) {
 r = i;
 for (int j = i+1; j < n; j++)
 if (fabs(a[j][i]) > fabs(a[r][i])) r = j;
 if (fabs(a[r][i]) < EPS) return false; // 为 0 无法化简
 if (r != i) swap(a[r], a[i]);
 for (int k = i+1; k < n; k++) {
 double f = a[k][i]/a[i][i];
 for (int j = i; j <= n; j++) a[k][j] -= f*a[i][j];
 }
 }

 // 回代
 for (int i = n-1; i >= 0; i--) {
 for (int j = i+1; j < n; j++)
 a[i][n] -= a[j][n]*a[i][j];
 a[i][n] /= a[i][i];
 }
 return true;
}
int main()
{
 cin >> n;
 for (int i = 0; i < n; i++)
 for (int j = 0; j <= n; j++)
 cin >> a[i][j];
 if (gauss()) {
 for (int i = 0; i < n; i++)
 printf("%.2lf\n", a[i][n]);
 } else {
 puts("No Solution");
 }
 return 0;
}

高斯-约旦消元

高斯-乔丹？这个消元法不需要回代，代码更简单。

将高斯消元的上三角形式化为对角线形式，就是高斯-约旦消元。最终结果的系数部分的矩阵只有对角线上的可以为 $1$，其余都为 $0$。在代码上的区别就是将回代的部分转移到前面的循环而已。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

bool gauss_jordan() {
 int r;
 for (int i = 0; i < n; i++) {
 r = i;
 for (int j = i+1; j < n; j++)
 if (fabs(a[j][i]) > fabs(a[r][i])) r = j;
 if (fabs(a[r][i]) < EPS) return false; // 为 0 无法化简
 if (r != i) swap(a[r], a[i]);
 for (int k = 0; k < n; k++) { // 区别：高斯消元不处理小于 i 的行，因此最后需要回代
 if (k == i) continue;
 double f = a[k][i]/a[i][i];
 for (int j = i; j <= n; j++) a[k][j] -= f*a[i][j];
 }
 }
 return true;
}
// 最后的输出部分改为：
printf("%.2lf\n", a[i][n]/a[i][i]);