有負權圖上的最短路算法 (Goldberg, 1995)

最近聽說有了一個有負權圖上的 \(O(m\log^8 m \log w)\) 算法, 感覺非常厲害, 於是覺得先來研讀一個早些的工作. 如果有可能的話再嘗試研讀新的算法!

注意: 本人主要是向論文學習了關鍵性質, 細節部分可能有所不同.

我們知道, OI 中常用的在負權圖上的 Bellman–Ford 算法可以在 \(O(nm)\) 時間內計算一個有負權圖的單源最短路徑, 或者確定這張圖存在一個負環. 而這篇文章給出了一個 \(O(\sqrt n m \log w)\) 的做法, 其中 \(w\) 是負權邊的絕對值上界.

pricing

這是本文涉及的第一個技術, 其實也就是大家熟知的最短路重賦權方法. 如果我們有一個 \(p\) 數組, 滿足

\[p(v) \leq p(u) + w(u, v), \]

那么新邊權

\[w_p(u, v) = w(u, v) + p(u) - p(v) \]

是非負的, 我們可以使用經典的 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) = w(u, v) + 2(d(u) - d(v)), \]

則必有 \(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}\) 倍, 那么經過

\[\sum_{i\geq 0} \sqrt{\frac{|K|}{\mathrm{e}^i}} \leq \frac 1{1 - \mathrm {e}^{-1/2}} \sqrt{|K|} \]

輪之后一定完成.

等等, 我們是不是還得夠快的找到這么一條鏈/反鏈啊?

在 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;
}