有负权图上的最短路算法 (Goldberg, 1995)

最近听说有了一个有负权图上的 O(mlog8mlogw) 算法, 感觉非常厉害, 于是觉得先来研读一个早些的工作. 如果有可能的话再尝试研读新的算法!

注意: 本人主要是向论文学习了关键性质, 细节部分可能有所不同.

我们知道, OI 中常用的在负权图上的 Bellman–Ford 算法可以在 O(nm) 时间内计算一个有负权图的单源最短路径, 或者确定这张图存在一个负环. 而这篇文章给出了一个 O(nmlogw) 的做法, 其中 w负权边的绝对值上界.

pricing

这是本文涉及的第一个技术, 其实也就是大家熟知的最短路重赋权方法. 如果我们有一个 p 数组, 满足

p(v)p(u)+w(u,v),

那么新边权

wp(u,v)=w(u,v)+p(u)p(v)

是非负的, 我们可以使用经典的 Dijkstra 算法求解. 显然在没有负环的时候, 这样的 p 数组是存在的, 因为最短路数组就满足这个条件.

但是这一转化的另一点好处在于, 我们求出这一 p 数组的条件是宽了很多的, 具体的设计将在后面看到.

scaling

这是本文涉及的第二个技术, 这个对于学过弱多项式复杂度费用流的同学来说可能比较熟悉了. 对于一张图 G=(V,E) 上的权值 w:EZ, 显然将所有边权除以 2 (不考虑取整的时候) 不影响答案和负环的存在性, 然后将权值增加, 不存在负环的图肯定依然不存在负环. 那么我们考虑先对于权值为 w=w/2 的图计算出 p 数组, 首先如果 w 里面所有数都是偶数, 那就有 w=2w, 自然 2p 就满足要求, 但实际上有些要将 2w 减去 1, 我们就是需要对此进行解决. 我们知道, 根据 p 数组的定义, 我们必然有 d(v)w(u,v)+d(u), 也即 w(u,v)+d(u)d(v)0. 我们考虑重赋权

wp(u,v)=w(u,v)+2(d(u)d(v)),

则必有 wp(u,v)1. 那么我们的目标就是: 通过对 p 进行适当的调整, 让 wp0. 我们将问题变成了 logw 轮满足 w1 的情况. 那么只要我们设计出一个对此在 O(nm) 时间内调整 price 的方法, 就完成了最初的承诺.

简单算法

我们先把目前 w0 的边拿出来, 称其为 Gp. 若 Gp 的 SCC 里有负边, 那这其实就直接说明了负环的存在性, 我们可以直接报告了. 接下来, 我们考虑缩点后的图. 我们称一个顶点集合是的当每个 vS 都有 v 通过 Gp 可达的边都还在 S 内部, 我们观察如果将这样一个 S 中的所有点的 p 值减去 1, 会发生什么.

  • 对于 S 内部的边, 显然没有影响, 因为两端的 p 值变化量相同.
  • 对于 S 连出的边的边权会减少 1, 但是注意到连出的边此时皆为正的, 这不会产生新的 1 边, 但有可能会加入新的 0 边.
  • 对于连入 S 的边, 边权加 1, 这说明连入 S 的负边会全部消失, 这是这个操作的关键用途.

如果我们将 S 取为某个点 v 通过 Gp 能到达的所有点 (除去其自身), 那么操作之后, v 点连出的所有负边都会被解决, 我们只需要 O(n) 轮就可以解决所有的 1 边了!

不幸的是, 每次操作之后, 图的结构会有不小的改变, 因为正权边变小之后会产生一些新的 0 边, 我们需要继续缩点判负环. 这样的话我们就得到了一个 O(nm) 的单轮算法.

但事实上, 我们距离最终的算法已相去不远, 这只需要我们用到一点组合观察, 并针对性的设计上述的操作方式.

最终算法

我们记有 1 作为出边的点集为 K, 在 Gp 上的"可达性"作为偏序, 那么我们考虑能否通过一种精心设计的操作一次性干掉 K 中的多个点.

注意: 为了后续论述的正确性, 这里的可达性需要定义为: u 可以经过一条 1 边, 然后经过 Gp 中的点到达 v.

反链

对于 K 中的一条反链 (antichain), 也即这些点两两不能到达, 那么令 S 为这些点在 V 上到达的点集之并, 我们对 S 进行操作, 那么这些点的出边都会被得到解决. 这和上述提到的做法并无不同, 复杂度为 O(m).

对于 K 中的一条链 (chain), 也即顶点 v1,v2,,vr 满足 vi 可以通过 Gp 到达 vi+1, 记 Sv 是 "v 直接通过 1 边到达的点在 Gp 中的可达点集之并", 那么我们考虑如果按照 i 从大到小的顺序操作, 会发生什么.

在还完全没有进行操作的时候, 显然我们是有 Sv1Sv2Svr 的, 那么一开始我们对 Svr 进行操作, 这会导致有些边的边权发生变化, 有的由正变 0, 也有的由 0 变正. 更靠前的 Svi 有可能变大的同时变小, 但它依然包含现在的 Svi+1, 这是因为它依旧能到达 vi+1, 并且原先 vi+1 连出的 1 边现在是 0, 仍在 Gp 内部.

但是有没有可能我们操作 S 之后 vi 并没被解决呢? 注意, 这个情况说明 Svi+1 里面有点可以到达 vi, 那负环就出现了.

这个操作相对而言维护起来就需要一点技巧了, 需要动态检查哪些边降为了 0, 这可以通过基数排序做到 (这些边一开始的权值不能超过 n, 所以基数排序和值域无关). 这样就做到了 O(m).

Dilworth

现在你有一个链上的算法, 你有一个反链上的算法. 现在让我们来过最后一关.

由 Dilworth 定理, 链和反链必有一者 |K|. 因此迭代 O(|K|) 轮后我们就完成了.

什么, 你问轮数为啥是 O(|K|)?

经过 |K| 轮, 集合大小每轮乘以一个不超过 11|K| 的因子, 那么集合变成原来的 1/e 倍, 那么经过

i0|K|ei11e1/2|K|

轮之后一定完成.

等等, 我们是不是还得够快的找到这么一条链/反链啊?

在 DAG 上 DP 可以求出每个点连出的最长链, 然后对所有极大点构成的反链进行检查即可, 这样即可在 O(m) 时间内找到.

综上, 一轮总共时间为 O(nm), 算法的总时间为 O(nmlogw).

实现

如下代码可以通过 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;
}
posted @   EntropyIncreaser  阅读(3450)  评论(12编辑  收藏  举报
相关博文:
阅读排行:
· 如何给本地部署的DeepSeek投喂数据,让他更懂你
· 超详细,DeepSeek 接入PyCharm实现AI编程!(支持本地部署DeepSeek及官方Dee
· 用 DeepSeek 给对象做个网站,她一定感动坏了
· .NET 8.0 + Linux 香橙派,实现高效的 IoT 数据采集与控制解决方案
· DeepSeek处理自有业务的案例:让AI给你写一份小众编辑器(EverEdit)的语法着色文件
点击右上角即可分享
微信分享提示