__gnu_cxx::rope的真面目——浅谈可持久化并查集测试

发表于 更新于 阅读次数： 本文字数： 11k 阅读时长 ≈ 10 分钟

最近又遇到了可以用可持久化并查集解决的问题，尝试用__gnu_cxx::rope水过，结果被卡成暴力分。鉴于已经有前车之鉴——一次模拟赛也是这样，但是BZOJ上的可持久化并查集却能水过，于是我打算做一些实验。

可持久化数组

由于可持久化并查集的基础是可持久化数组，因此我设计了一道测试用的模板题。

为了简单起见，给定数组$a[1\dots n]$，进行$m$次修改操作，然后$q$个询问，求某次修改后某个位置上的数。由于仅用于测试速度，所以数据随机。$0\le a[i]\le 10^9$。

1
2
3
4
5
6
7
8
9
10
11
12
13
5
1 2 3 4 5
4
2 6
3 2
2 3
2 1
5
1 3
4 3
3 5
4 2
2 2
1
2
3
4
5
3
2
5
1
6
id $n,m,q$
1 $10^4$
2 $3\times 10^4$
3 $10^5$
4 $3\times 10^5$
5 $10^6$
gen.cpp
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
#include <fstream>
#include <random>
#include <ctime>
using namespace std;
ofstream fout("rope.in");
const int n = 3e5, m = 3e5, q = 3e5, val = 1e9;
minstd_rand gen(time(NULL));
int random(int l, int r)
{
	return uniform_int_distribution<>(l, r)(gen);
}
int main()
{
	fout << n << endl;
	for (int i = 1; i <= n; i++)
		fout << random(0, val) << ' ';
	fout << endl
		 << m << endl;
	for (int i = 1; i <= m; i++)
		fout << random(1, n) << ' ' << random(0, val) << endl;
	fout << q << endl;
	for (int i = 1; i <= q; i++)
		fout << random(1, m) << ' ' << random(1, n) << endl;
	return 0;
}

测试在Manjaro下进行，不计输入输出时间，估计折算净内存使用（有较大的误差，仅作参考）。即用CCR-Plus的monitor来获取内存使用（iorun不支持），因此设置限制为“10 s, 2 GB”。在使用前，需要把monitor代码中的“进程被阻塞”判定删掉。

__gnu_cxx::rope

定义在<ext/rope>中。可以看出，这和pb_ds一样是一个GNU库，因此只有libstdc++中有。据说也是由SGI开发的，SGI更为著名的是STL。这个模板类原本是为了快速字符串编辑，像是不少古老的文本编辑器题目可以用rope水过。但也可以把存储类型改成int，再加上特殊的存储方式，使复制的复杂度低于线性，于是就成了“可持久化数组”。

使用方法可以查阅man __gnu_cxx::rope，这可能是用NOI Linux的优势，可以随便查函数。为了实现快速复制，一般用指针管理更快，调用复制构造就可以快速复制了。rope的operator []和at()都不能作为左值，修改一个数要用replace()方法。

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 <fstream>
#include <ext/rope>
#include <ctime>
#include <iostream>
using namespace std;
using __gnu_cxx::rope;
ifstream fin("rope.in");
ofstream fout("rope.out");
const int N = 1000005, M = 1000005, Q = 1000005;
int a[N], ans[Q];
pair<int, int> mod[M], que[Q];
rope<int> *root[M];
int main()
{
	int n;
	fin >> n;
	for (int i = 1; i <= n; i++)
		fin >> a[i];
	int m;
	fin >> m;
	for (int i = 1; i <= m; i++)
		fin >> mod[i].first >> mod[i].second;
	int q;
	fin >> q;
	for (int i = 1; i <= q; i++)
		fin >> que[i].first >> que[i].second;
	clock_t t = clock();
	root[0] = new rope<int>(a, a + n + 1);//比一个一个加入更高效
	for (int i = 1; i <= m; i++)
	{
		root[i] = new rope<int>(*root[i - 1]);
		root[i]->replace(mod[i].first, mod[i].second);
	}
	for (int i = 1; i <= q; i++)
		ans[i] = root[que[i].first]->at(que[i].second);
	cerr << clock() - t << endl;
	for (int i = 1; i <= q; i++)
		fout << ans[i] << endl;
	return 0;
}
id O0 time O2 time mem
1 81 ms 43 ms 30.6 MB
2 302 ms 157 ms 107.37 MB
3 1184 ms 639 ms 425.05 MB
4 4313 ms 2227 ms 1425.61 MB
5 N/A N/A N/A

librope

宣传的比SGI rope更高效的实现，参见https://gist.github.com/josephg/3474848#file-results。不过仔细观察就可以发现，librope是用跳表实现的，而不是树，而且测试的项目是大量插入和少量删除操作。

不幸的是，librope是用C写的，不支持模板，只支持字符。更讨厌的是只能用C风格字符串初始化和修改，我不得不把int转为4个字符，而且还得保证中间没有0，一种简单的方法是用255进制。

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
#include <fstream>
#include "rope.h"
#include <ctime>
#include <iostream>
#include <cstdint>
using namespace std;
ifstream fin("rope.in");
ofstream fout("rope.out");
const int N = 1000005, M = 1000005, Q = 1000005;
int a[N], ans[Q];
pair<int, int> mod[M], que[Q];
rope *root[M];
int main()
{
	int n;
	fin >> n;
	for (int i = 1; i <= n; i++)
		fin >> a[i];
	int m;
	fin >> m;
	for (int i = 1; i <= m; i++)
		fin >> mod[i].first >> mod[i].second;
	int q;
	fin >> q;
	for (int i = 1; i <= q; i++)
		fin >> que[i].first >> que[i].second;
	clock_t t = clock();
	uint8_t *as = new uint8_t[n * 4 + 1];
	for (int i = 1; i <= n; i++)
		for (int j = 0; j < 4; j++)
		{
			as[(i - 1) * 4 + j] = a[i] % 255 + 1;
			a[i] /= 255;
		}
	as[n * 4] = 0;
	root[0] = rope_new_with_utf8(as);
	for (int i = 1; i <= m; i++)
	{
		root[i] = rope_copy(root[i - 1]);
		uint8_t as[5];
		for (int j = 0; j < 4; j++)
		{
			as[j] = mod[i].second % 255 + 1;
			mod[i].second /= 255;
		}
		as[4] = 0;
		rope_del(root[i], (mod[i].first - 1) * 4, 4);
		rope_insert(root[i], (mod[i].first - 1) * 4, as);
	}
	for (int i = 1; i <= q; i++)
	{
		rope *tmp = rope_copy(root[que[i].first]);
		rope_del(tmp, 0, (que[i].second - 1) * 4);
		rope_del(tmp, 4, (n - que[i].second) * 4);
		uint8_t *as = rope_create_cstr(tmp);
		for (int j = 3; j >= 0; j--)
			ans[i] = ans[i] * 255 + as[j] - 1;
	}
	cerr << clock() - t << endl;
	for (int i = 1; i <= q; i++)
		fout << ans[i] << endl;
	return 0;
}

注意要先用gcc编译librope，再用g++编译这份代码。测试结果悲惨的证明，根本不支持快速复制，看来上当了！

id O0 time O2 time mem
1 652 ms 574 ms 680.71 MB
2 N/A N/A N/A

主席树

最为简单的可持久化数据结构，用作可持久化数组时，非叶节点不用保存值。其实可以把叶节点的值放在左儿子（或右儿子）来节省不少空间。空间复杂度和倍增一样达到$\mathcal O(n\log n)$，每次操作都是$\mathcal O(\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
85
#include <fstream>
#include <ctime>
#include <iostream>
using namespace std;
ifstream fin("rope.in");
ofstream fout("rope.out");
const int N = 1000005, M = 1000005, Q = 1000005;
int a[N], ans[Q], cnt;
pair<int, int> mod[M], que[Q];
typedef int ptr;//容易区分是下标还是值
struct node
{
	int ls, rs;//这个没办法，相当于手动union
} tree[M * 25];
ptr root[M];
void build(ptr &id, int l, int r)
{
	id = ++cnt;
	if (l == r)
		tree[id].ls = a[l];
	else
	{
		int mid = (l + r) / 2;
		build(tree[id].ls, l, mid);
		build(tree[id].rs, mid + 1, r);
	}
}
void modify(ptr &id, ptr ref, int l, int r, int x, int val)
{
	id = ++cnt;
	if (l == r)
		tree[id].ls = val;
	else
	{
		int mid = (l + r) / 2;
		if (x <= mid)
		{
			modify(tree[id].ls, tree[ref].ls, l, mid, x, val);
			tree[id].rs = tree[ref].rs;//不修改的直接复制
		}
		else
		{
			modify(tree[id].rs, tree[ref].rs, mid + 1, r, x, val);
			tree[id].ls = tree[ref].ls;
		}
	}
}
int query(ptr id, int l, int r, int x)
{
	if (l == r)
		return tree[id].ls;
	int mid = (l + r) / 2;
	if (x <= mid)
		return query(tree[id].ls, l, mid, x);
	return query(tree[id].rs, mid + 1, r, x);
}
int main()
{
	int n;
	fin >> n;
	for (int i = 1; i <= n; i++)
		fin >> a[i];
	int m;
	fin >> m;
	for (int i = 1; i <= m; i++)
		fin >> mod[i].first >> mod[i].second;
	int q;
	fin >> q;
	for (int i = 1; i <= q; i++)
		fin >> que[i].first >> que[i].second;
	clock_t t = clock();
	build(root[0], 1, n);
	for (int i = 1; i <= m; i++)
	{
		root[i] = root[i - 1];
		modify(root[i], root[i], 1, n, mod[i].first, mod[i].second);
	}
	for (int i = 1; i <= q; i++)
		ans[i] = query(root[que[i].first], 1, n, que[i].second);
	cerr << clock() - t << endl;
	//cerr << cnt << endl;
	for (int i = 1; i <= q; i++)
		fout << ans[i] << endl;
	return 0;
}
id O0 time O2 time
1 5 ms 7 ms
2 18 ms 18 ms
3 92 ms 56 ms
4 334 ms 215 ms
5 1308 ms 871 ms

由于开的是静态数组，内存估计恒定为195 MB。

分块

显然可以做，但是好像复杂度不平衡。

留坑待填

Random Access List

参见https://www.cs.oberlin.edu/~jwalker/ra-list/，一种没听说过的数据结构，还可以做我以前出的数列。

留坑待填

统计

可持久化并查集

按秩合并

一种常见的伎俩，即合并时把少的合并到多的，常称为启发式合并。可以保证并查集的树高为$\mathcal O(\log n)$级别。若还是$m$次操作$q$次询问，复杂度$\mathcal O((m+q)\log^2 n)$。

路径压缩？

据说不能路径压缩，但是其实好像也没有坏处，但会慢不少。在按秩合并的基础上，复杂度$\mathcal O(\alpha m\log n+q\log^2 n)$，因为询问的高度最坏可以达到$\mathcal O(\log n)$。

例题

不要认为这个没什么用，其实有时候可以暴力过一些题。

模板

如果有空，可以用上面的方法去试一下：洛谷P3919、洛谷P3402，当然还有著名的BZOJ 3673、3674。

NOIp2013 货车运输

大意是给定一个无向图，多次询问两点之间的路径，最大化路径上最小边。标准做法是建最大生成树，然后树上倍增。最近才做了这题，然后就开始研究可持久化并查集。

但是我一开始没有想到只有最大生成树上的边有效，于是开始折腾。首先注意到可以二分答案，在并查集上连上大于等于答案的边，就可以判断连通性。复杂度$\mathcal O(qm\log z)$。然后发现不用二分，离线处理，直接从大到小加边，每次判断所有询问。复杂度$\mathcal O(qm)$，这样就是暴力的60分。

还能怎么改进呢？如果用可持久化并查集，就能保存每次加边后的并查集，询问只要二分判断连通即可。然后我就用了SGI rope，结果还是只有60分，剩下的MLE。

我不肯就此罢休，于是改成了主席树。不幸的是还是60分，剩下的RE，恰好是满数据。主席树的节点应该够了呀。搞来大数据才发现其实是写错了，但只有询问足够大才能遇到这种错误。还好这种错误我在写的时候就注意到了，由于我用了路径压缩，如果询问直接调用getf()就可能会修改前面的版本，而造成后面的版本不兼容。因此要么先复制一份，要么写一个“只读”版本的getf()。加上按秩合并优化常数，最后终于过了。感觉完全是在作死，比赛不太可能写出来。

复杂度$\mathcal O(\alpha m\log n+q\log^3 n)$，好像很慢的样子，其实后面的3次log不会满的，大约只有2次。但在当年可能还是会被卡常，毕竟标算只有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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
#include <iostream>
#include <algorithm>
using namespace std;
const int N = 10005, M = 50005, NODES = 5e6;
struct edge
{
	int u, v, w;
	bool operator<(const edge &rhs) const
	{
		return w > rhs.w;
	}
} e[M];
typedef int ptr;
struct node
{
	int ls, rs;
} tree[NODES];
ptr root[M];
int n, m, cnt, sz[N];
void build(ptr &id, int l, int r)
{
	id = ++cnt;
	if (l == r)
		tree[id].ls = l;
	else
	{
		int mid = (l + r) / 2;
		build(tree[id].ls, l, mid);
		build(tree[id].rs, mid + 1, r);
	}
}
void modify(ptr &id, ptr rt, int l, int r, int x, int val)
{
	id = ++cnt;
	if (l == r)
		tree[id].ls = val;
	else
	{
		int mid = (l + r) / 2;
		if (x <= mid)
		{
			modify(tree[id].ls, tree[rt].ls, l, mid, x, val);
			tree[id].rs = tree[rt].rs;
		}
		else
		{
			modify(tree[id].rs, tree[rt].rs, mid + 1, r, x, val);
			tree[id].ls = tree[rt].ls;
		}
	}
}
int query(ptr id, int l, int r, int x)
{
	if (l == r)
		return tree[id].ls;
	int mid = (l + r) / 2;
	if (x <= mid)
		return query(tree[id].ls, l, mid, x);
	return query(tree[id].rs, mid + 1, r, x);
}
int getf(ptr &id, int x)
{
	int fat = query(id, 1, n, x);
	if (fat == x)
		return x;
	int rt = getf(id, fat);
	modify(id, id, 1, n, x, rt);
	return rt;
}
int getf_readonly(ptr id, int x)
{
	int fat = query(id, 1, n, x);
	if (fat == x)
		return x;
	return getf_readonly(id, fat);
}
int main()
{
	ios::sync_with_stdio(false);
	cin >> n >> m;
	for (int i = 1; i <= m; i++)
		cin >> e[i].u >> e[i].v >> e[i].w;
	sort(e + 1, e + m + 1);
	build(root[0], 1, n);
	fill(sz + 1, sz + n + 1, 1);
	for (int i = 1; i <= m; i++)
	{
		root[i] = root[i - 1];
		int ru = getf(root[i], e[i].u), rv = getf(root[i], e[i].v);//改成getf_readonly更快
		if (ru != rv)
		{
			if (sz[ru] > sz[rv])
			{
				modify(root[i], root[i], 1, n, rv, ru);
				sz[ru] += sz[rv];
			}
			else
			{
				modify(root[i], root[i], 1, n, ru, rv);
				sz[rv] += sz[ru];
			}
		}
	}
	int q;
	cin >> q;
	while (q--)
	{
		int u, v;
		cin >> u >> v;
		int l = 1, r = m;
		while (l < r)
		{
			int mid = (l + r) / 2;
			if (getf_readonly(root[mid], u) == getf_readonly(root[mid], v))
				r = mid;
			else
				l = mid + 1;
		}
		if (getf_readonly(root[l], u) != getf_readonly(root[l], v))
			cout << -1 << endl;
		else
			cout << e[l].w << endl;
	}
	return 0;
}

SillyHook’s history

随便取了个名字，我们这届的2017-07-17模拟赛T3（去年暑假集训），也有公开的https://wenku.baidu.com/view/16dc09b8250c844769eae009581b6bd97e19bc74.html。

当年直接用了SGI rope，结果卡成30分，本来可以写写其他部分分的。