有负权图上的最短路算法 (Goldberg, 1995)
最近听说有了一个有负权图上的 \(O(m\log^8 m \log w)\) 算法, 感觉非常厉害, 于是觉得先来研读一个早些的工作. 如果有可能的话再尝试研读新的算法!
注意: 本人主要是向论文学习了关键性质, 细节部分可能有所不同.
我们知道, OI 中常用的在负权图上的 Bellman–Ford 算法可以在 \(O(nm)\) 时间内计算一个有负权图的单源最短路径, 或者确定这张图存在一个负环. 而这篇文章给出了一个 \(O(\sqrt n m \log w)\) 的做法, 其中 \(w\) 是负权边的绝对值上界.
pricing
这是本文涉及的第一个技术, 其实也就是大家熟知的最短路重赋权方法. 如果我们有一个 \(p\) 数组, 满足
那么新边权
是非负的, 我们可以使用经典的 Dijkstra 算法求解. 显然在没有负环的时候, 这样的 \(p\) 数组是存在的, 因为最短路数组就满足这个条件.
但是这一转化的另一点好处在于, 我们求出这一 \(p\) 数组的条件是宽了很多的, 具体的设计将在后面看到.
scaling
这是本文涉及的第二个技术, 这个对于学过弱多项式复杂度费用流的同学来说可能比较熟悉了. 对于一张图 \(G=(V,E)\) 上的权值 \(w: E \to \mathbb Z\), 显然将所有边权除以 \(2\) (不考虑取整的时候) 不影响答案和负环的存在性, 然后将权值增加, 不存在负环的图肯定依然不存在负环. 那么我们考虑先对于权值为 \(w' = \lceil w / 2 \rceil\) 的图计算出 \(p\) 数组, 首先如果 \(w\) 里面所有数都是偶数, 那就有 \(w = 2w'\), 自然 \(2p\) 就满足要求, 但实际上有些要将 \(2w'\) 减去 \(1\), 我们就是需要对此进行解决. 我们知道, 根据 \(p\) 数组的定义, 我们必然有 \(d(v) \leq w'(u, v) + d(u)\), 也即 \(w'(u,v) + d(u) - d(v) \geq 0\). 我们考虑重赋权
则必有 \(w_p(u, v) \geq -1\). 那么我们的目标就是: 通过对 \(p\) 进行适当的调整, 让 \(w_p\geq 0\). 我们将问题变成了 \(\log w\) 轮满足 \(w\geq -1\) 的情况. 那么只要我们设计出一个对此在 \(O(\sqrt n m)\) 时间内调整 price 的方法, 就完成了最初的承诺.
简单算法
我们先把目前 \(w \leq 0\) 的边拿出来, 称其为 \(G_p\). 若 \(G_p\) 的 SCC 里有负边, 那这其实就直接说明了负环的存在性, 我们可以直接报告了. 接下来, 我们考虑缩点后的图. 我们称一个顶点集合是闭的当每个 \(v\in S\) 都有 \(v\) 通过 \(G_p\) 可达的边都还在 \(S\) 内部, 我们观察如果将这样一个 \(S\) 中的所有点的 \(p\) 值减去 \(1\), 会发生什么.
- 对于 \(S\) 内部的边, 显然没有影响, 因为两端的 \(p\) 值变化量相同.
- 对于 \(S\) 连出的边的边权会减少 \(1\), 但是注意到连出的边此时皆为正的, 这不会产生新的 \(-1\) 边, 但有可能会加入新的 \(0\) 边.
- 对于连入 \(S\) 的边, 边权加 \(1\), 这说明连入 \(S\) 的负边会全部消失, 这是这个操作的关键用途.
如果我们将 \(S\) 取为某个点 \(v\) 通过 \(G_p\) 能到达的所有点 (除去其自身), 那么操作之后, \(v\) 点连出的所有负边都会被解决, 我们只需要 \(O(n)\) 轮就可以解决所有的 \(-1\) 边了!
不幸的是, 每次操作之后, 图的结构会有不小的改变, 因为正权边变小之后会产生一些新的 \(0\) 边, 我们需要继续缩点判负环. 这样的话我们就得到了一个 \(O(nm)\) 的单轮算法.
但事实上, 我们距离最终的算法已相去不远, 这只需要我们用到一点组合观察, 并针对性的设计上述的操作方式.
最终算法
我们记有 \(-1\) 作为出边的点集为 \(K\), 在 \(G_p\) 上的"可达性"作为偏序, 那么我们考虑能否通过一种精心设计的操作一次性干掉 \(K\) 中的多个点.
注意: 为了后续论述的正确性, 这里的可达性需要定义为: \(u\) 可以经过一条 \(-1\) 边, 然后经过 \(G_p\) 中的点到达 \(v\).
反链
对于 \(K\) 中的一条反链 (antichain), 也即这些点两两不能到达, 那么令 \(S\) 为这些点在 \(V\) 上到达的点集之并, 我们对 \(S\) 进行操作, 那么这些点的出边都会被得到解决. 这和上述提到的做法并无不同, 复杂度为 \(O(m)\).
链
对于 \(K\) 中的一条链 (chain), 也即顶点 \(v_1, v_2, \dots, v_r\) 满足 \(v_i\) 可以通过 \(G_p\) 到达 \(v_{i+1}\), 记 \(S_v\) 是 "\(v\) 直接通过 \(-1\) 边到达的点在 \(G_p\) 中的可达点集之并", 那么我们考虑如果按照 \(i\) 从大到小的顺序操作, 会发生什么.
在还完全没有进行操作的时候, 显然我们是有 \(S_{v_1} \supset S_{v_2} \supset \dots \supset S_{v_r}\) 的, 那么一开始我们对 \(S_{v_r}\) 进行操作, 这会导致有些边的边权发生变化, 有的由正变 \(0\), 也有的由 \(0\) 变正. 更靠前的 \(S_{v_i}\) 有可能变大的同时变小, 但它依然包含现在的 \(S_{v_{i+1}}\), 这是因为它依旧能到达 \(v_{i+1}\), 并且原先 \(v_{i+1}\) 连出的 \(-1\) 边现在是 \(0\), 仍在 \(G_p\) 内部.
但是有没有可能我们操作 \(S\) 之后 \(v_i\) 并没被解决呢? 注意, 这个情况说明 \(S_{v_{i+1}}\) 里面有点可以到达 \(v_i\), 那负环就出现了.
这个操作相对而言维护起来就需要一点技巧了, 需要动态检查哪些边降为了 \(0\), 这可以通过基数排序做到 (这些边一开始的权值不能超过 \(n\), 所以基数排序和值域无关). 这样就做到了 \(O(m)\).
Dilworth
现在你有一个链上的算法, 你有一个反链上的算法. 现在让我们来过最后一关.
由 Dilworth 定理, 链和反链必有一者 \(\geq \sqrt{|K|}\). 因此迭代 \(O(\sqrt {|K|})\) 轮后我们就完成了.
什么, 你问轮数为啥是 $O(\sqrt{|K|})$?
经过 \(\sqrt{|K|}\) 轮, 集合大小每轮乘以一个不超过 \(1 - \frac 1{\sqrt{|K|}}\) 的因子, 那么集合变成原来的 \(\leq 1/\mathrm{e}\) 倍, 那么经过
轮之后一定完成.
等等, 我们是不是还得够快的找到这么一条链/反链啊?
在 DAG 上 DP 可以求出每个点连出的最长链, 然后对所有极大点构成的反链进行检查即可, 这样即可在 \(O(m)\) 时间内找到.
综上, 一轮总共时间为 \(O(\sqrt n m)\), 算法的总时间为 \(O(\sqrt n m \log w)\).
实现
如下代码可以通过 LuoguP5960 差分约束算法.
#include <bits/stdc++.h>
using namespace std;
const int V = 5005, E = 5005;
int N, M;
int u[E], v[E], w[E];
int p[V];
bool mark[V];
vector<pair<int, int>> G[V], G2[V];
int t, N2;
stack<int> stk;
int dfn[V], low[V], belongs[V], topo[V];
void tarjan(int u) {
dfn[u] = low[u] = ++t; stk.push(u);
for (const auto &e : G[u]) {
int v = e.first, w = e.second + p[u] - p[v];
if (w > 0) continue;
if (!dfn[v]) {
tarjan(v);
low[u] = min(low[u], low[v]);
} else if (belongs[v] == -1) {
low[u] = min(low[u], dfn[v]);
}
}
if (dfn[u] == low[u]) {
topo[N2++] = u;
int v;
do {
v = stk.top(); stk.pop();
belongs[v] = u;
} while (v != u);
}
}
bool SCC() {
t = N2 = 0;
memset(dfn, 0, sizeof(dfn));
memset(mark, 0, sizeof(mark));
memset(belongs, -1, sizeof(belongs));
for (int i = 0; i != N; ++i) G2[i].clear();
for (int i = 0; i != N; ++i) if (!dfn[i]) tarjan(i);
#ifdef ELEGIA_DEBUG
int marks = 0;
#endif // ELEGIA_DEBUG
for (int i = 0; i != N; ++i) for (const auto &e : G[i]) {
int j = e.first, w = e.second + p[i] - p[j];
if (w < 0) {
if (belongs[i] == belongs[j]) return true;
mark[belongs[i]] = true;
#ifdef ELEGIA_DEBUG
++marks;
#endif // ELEGIA_DEBUG
}
G2[belongs[i]].emplace_back(belongs[j], w);
}
#ifdef ELEGIA_DEBUG
cerr << "marks: " << marks << endl;
#endif // ELEGIA_DEBUG
return false;
}
int vis[V], pref[V], sec[V], len[V], chn[V];
vector<int> delay[V];
void chain(int start) {
int len = 0;
while (start != -1) {
chn[++len] = start;
start = sec[start];
}
memset(vis, 0, sizeof(vis));
queue<int> q;
for (int i = len; i; --i) {
int s = chn[i];
for (const auto &e : G2[s]) {
int v = e.first, w = e.second;
if (vis[v] || w >= 0) continue;
vis[v] = i;
q.push(v);
}
for (int u : delay[i]) if (!vis[u]) {
vis[u] = i;
q.push(u);
}
delay[i].clear();
while (!q.empty()) {
int u = q.front(); q.pop();
for (const auto &e : G2[u]) {
int v = e.first, w = e.second;
if (w > 0) {
if (w < i) delay[i - w].push_back(v);
continue;
}
if (vis[v]) continue;
vis[v] = i;
q.push(v);
}
}
}
}
bool refine() {
while (true) {
{
bool flag = false;
for (int u = 0; u != N; ++u) {
for (const auto &e : G[u]) {
int v = e.first, w = e.second + p[u] - p[v];
if (w < 0) {
flag = true;
break;
}
}
if (flag) break;
}
if (!flag) return false;
}
if (SCC()) return true;
memset(len, 0, sizeof(len));
memset(pref, -1, sizeof(pref));
memset(sec, -1, sizeof(sec));
int longest = -1;
for (int i = 0; i != N2; ++i) {
int u = topo[i];
for (const auto &e : G2[u]) {
int v = e.first, w = e.second;
if (w > 0) continue;
if (pref[v] != -1 && len[v] > len[u]) {
len[u] = len[v];
pref[u] = pref[v];
}
}
if (mark[u]) {
int opt = 0;
for (const auto &e : G2[u]) {
int v = e.first, w = e.second;
if (w >= 0) continue;
if (pref[v] != -1 && len[v] > opt) {
opt = len[v];
sec[u] = pref[v];
}
}
if (++opt > len[u]) {
len[u] = opt;
pref[u] = u;
} else sec[u] = -1;
}
if (longest == -1 || len[u] > len[longest])
longest = u;
}
int size = 0;
queue<int> q;
memset(vis, 0, sizeof(vis));
for (int i = N2 - 1; i >= 0; --i) {
int u = topo[i];
if (vis[u] || !mark[u]) continue;
++size;
for (const auto &e : G2[u]) {
int v = e.first, w = e.second;
if (!vis[v] && w < 0) {
vis[v] = true;
q.push(v);
}
}
while (!q.empty()) {
int u = q.front(); q.pop();
for (const auto &e : G2[u]) {
int v = e.first, w = e.second;
if (!vis[v] && w <= 0) {
vis[v] = true;
q.push(v);
}
}
}
}
#ifdef ELEGIA_DEBUG
cerr << "antichain: " << size << ", chain: " << len[longest] << endl;
#endif // ELEGIA_DEBUG
if (size < len[longest]) chain(longest);
for (int i = 0; i != N; ++i)
p[i] -= vis[belongs[i]];
}
}
bool priceScaling() {
int L = 0;
for (int i = 0; i != M; ++i)
while ((1 << L) <= -w[i]) ++L;
memset(p, 0, sizeof(p));
for (int d = L - 1; d >= 0; --d) {
for (int i = 0; i != N; ++i) G[i].clear();
for (int i = 0; i != N; ++i) p[i] *= 2;
for (int i = 0; i != M; ++i)
G[u[i]].emplace_back(v[i], -((-w[i]) >> d));
if (refine()) return true;
#ifdef ELEGIA_DEBUG
cerr << "finish scaling 2^" << d << endl;
#endif // ELEGIA_DEBUG
}
return false;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr); cout.tie(nullptr);
cin >> N >> M;
for (int i = 0; i != M; ++i) {
cin >> v[i] >> u[i] >> w[i];
--u[i]; --v[i];
}
if (priceScaling()) {
cout << "NO\n";
} else for (int i = 0; i != N; ++i)
cout << p[i] << " \n"[i == N - 1];
return 0;
}
