莫队算法 (Mo's Algorithm)

省赛交了不熟莫队的学费之后，决定写篇博客复习一下。由于本人非常鄙视此类暴力算法(因为涉及分块，感觉很不优美，而且我分块姿势也不熟练)，于是一直没有重视，省赛就被教育了……

比如GDCPC2019广东省赛就有这么一道题：

　　给定n,m,k，一个长度为n的数组，m次询问。每次询问给出区间[l,r]，要求计算区间中有多少个a[i]、a[j]满足i<j && abs(a[i]-a[j])<=k。n,m<=27000, k<=1e9。TL=1s，OL=256mb。

当时第一反应是线段树，然而马上感觉非常不可做，因为线段树只能维护符合线性性的属性(比如区间最值、区间和之类的)；而这道题如果用线段树维护，等于要把左右子区间的信息再一次合并，左右子区间的信息完全没有贡献(因为不符合线性性，不能通过线性操作来维护当前节点的信息)，光是maintain的时间复杂度就已经达到了O(n^2*logn)，更不要说建树和查询了，完全fake做法。

那应该怎么做呢？由于这道题没有修改操作，所以可以考虑离线处理，最明显的处理区间问题的离线算法就是莫队了……莫队本质上是一种分块算法，通常用来离线解决只查询不修改的区间问题。

让我们看一道更水的题来学习莫队算法：

　　【例题】给定一个大小为N的数组，数组中所有元素的大小a[i]<=N。你需要回答M个查询。每个查询的形式是L，R。你需要回答在范围[ L,R ]这个区间内数字的出现次数刚好是k的数字种数。k<=n,m<=30000，1<=L<r<=n。TL=1s，OL=256mb。

这个题能用线段树做吗？依然不行。原因还是我上面提到的问题：在维护segT[currentPosition]的信息时，不能简单地通过segT[leftSon]和segT[rightSon]的信息来计算得出。每次维护segT[currentPosition]的信息时，都要重新处理左右子区间的信息来确定相互之间的影响，等于一棵fake线段树。

既然不能线段树，我们怎么思考这个问题呢？不妨假设我已经处理完了某个区间[l,r]的信息，现在我要处理区间[l,r+1]的信息，那我们就可以在O(1)时间内处理完毕：只需开一个cnt数组来记录数字的出现次数，再cnt [ a [ r+1 ] ] ++即可。这样，我们就可以知道区间[l±1,r±1]的信息，就可以用莫队算法了。

莫队算法的一种实现方式是：离线得到了一堆需要处理的区间后，合理地安排这些区间计算的次序以得到一个较优的复杂度。假设我们当前已知区间[l1,r1]的信息，要转移到区间[l2,r2]，由于l、r只能单步转移，那么时间复杂度为O(abs(l1-l2)+abs(r1-r2))。要是把这两个区间看作是平面上的两个整点，就变成了曼哈顿距离。整体的时间复杂度必然为整棵树的曼哈顿距离之和。显然当整棵树为MST时，时间复杂度最优。那么在求解答案时只要根据曼哈顿MST的dfs序求解即可。

还有另一种暴力写法：先对序列分块，然后以询问左端点所在的分块的序号为第一关键字，右端点的大小为第二关键字进行排序，按照排序好的顺序计算，复杂度就会大大降低。

分块相同时，右端点递增是 $O(N)$ 的，分块共有 $O(\sqrt{N} )$ 个，复杂度为 $O(N^{1.5} )$
分块转移时，右端点最多变化 $N$ ，分块共有 $O(\sqrt{N} )$ 个，复杂度为 $O(N^{1.5} )$
分块相同时，左端点最多变化 $\sqrt{N}$ ，分块转移时，左端点最多变化 $2\sqrt{N}$ ，共有 $N$ 个询问，复杂度为 $O(N^{1.5} )$

故总时间复杂度 $O(N^{1.5} )$

以例题为例。不妨给出一组样例：n=9 (数组内容忽略，不重要) ，m=8。查询区间为：[2,3], [1,4], [4,5], [1,6], [7,9], [8,9], [5,8], [6,8]。

因为最大范围不超过n，所以我们以(int)sqrt(n)为大小，对区间进行分块。在每一个块中，我们按r从小到大的顺序排列。所以上面的排序结果是：

{ (2,3) (1,4) (1,6)，(4,5) (5,8) (6,8)，(7,9) (8,9) }

这一步的排序可以这样实现：

1 unit=(int)sqrt(n);
2 bool operator<(const Interval &rhs)const
3 {
4     if (l / unit != rhs.l / unit) return l / unit < rhs.l / unit;
5     else return r < rhs.r;
6 }

考虑到在同一个块的时候，由于L的范围确定，故每次L的偏移量是O(sqrt(n))；但是R的范围没有确定，故R的偏移量是O(n)。

那么从一个块到另一个块呢? 显然我们不用考虑R的偏移量，依然是O(n)，而L明显最多也是2*sqrt(n)。在这种情况下，很快就会到下下一块。所以也是O(sqrt(n))。

由于有sqrt(n)个块，所以R的总偏移量是O(n*sqrt(n))，而M个询问，每个询问都可以让L偏移O(sqrt(n))，所以L的总偏移量O(m*sqrt(n))。

注意，R的偏移量和询问数目没有直接关系。而L则恰恰相反；L的偏移量我们刚才也说明了，和块的个数没有直接关系。所以总的时间复杂度是：O((n+m)*sqrt(n))。

还有道莫队经典题也分享一下：bzoj2038

大意是询问区间内任意选两个数为同一个数的概率并化为最简分数。

设在某一区间内共有颜色 $a1,a2,a3...an$ ，每双袜子的个数为 $b1,b2,b3...bn$ 。

答案为 $(\sum_{i=1}^{n}{b_{i}(b_{i}-1)/2} )/((R-L+1)(R-L)/2)$ 。

化简 $(\sum_{i=1}^{n}{b_{i}^{2} }-b)/((R-L+1)(R-L)/2)$ 。

即 $((\sum_{i=1}^{n}{b_{i}^{2} })-(R-L+1))/((R-L+1)(R-L)/2)$ 。

所以只需要用莫队处理每个区间内不同数字的平方和就好了。

分块写法：

 1 //分块
 2 #include <bits/stdc++.h>
 3 /* define */
 4 #define ll long long
 5 #define dou double
 6 #define pb emplace_back
 7 #define mp make_pair
 8 #define fir first
 9 #define sec second
10 #define sot(a,b) sort(a+1,a+1+b)
11 #define rep1(i,a,b) for(int i=a;i<=b;++i)
12 #define rep0(i,a,b) for(int i=a;i<b;++i)
13 #define repa(i,a) for(auto &i:a)
14 #define eps 1e-8
15 #define int_inf 0x3f3f3f3f
16 #define ll_inf 0x7f7f7f7f7f7f7f7f
17 #define lson curPos<<1
18 #define rson curPos<<1|1
19 /* namespace */
20 using namespace std;
21 /* header end */
22 
23 const int maxn = 5e4 + 10;
24 int n, m, unit, num[maxn], a[maxn];
25 struct Interval
26 {
27     int l, r, id;
28     bool operator<(const Interval &rhs)const
29     {
30         if (l / unit != rhs.l / unit) return l / unit < rhs.l / unit;
31         else return r < rhs.r;
32     }
33 } interval[maxn];
34 struct Ans
35 {
36     ll a, b;
37     void reduce()
38     {
39         ll g = __gcd(a, b);
40         a /= g, b /= g;
41     }
42 } ans[maxn];
43 
44 void solve()
45 {
46     ll tmp = 0;
47     rep0(i, 0, maxn) num[i] = 0;
48     int l = 1, r = 0; //初始区间
49     rep1(i, 1, m)
50     {
51         while (r < interval[i].r)
52         {
53             r++;
54             tmp -= (ll)num[a[r]] * num[a[r]];
55             num[a[r]]++;
56             tmp += (ll)num[a[r]] * num[a[r]];
57         }
58         while (r > interval[i].r)
59         {
60             tmp -= (ll)num[a[r]] * num[a[r]];
61             num[a[r]]--;
62             tmp += (ll)num[a[r]] * num[a[r]];
63             r--;
64         }
65         while (l < interval[i].l)
66         {
67             tmp -= (ll)num[a[l]] * num[a[l]];
68             num[a[l]]--;
69             tmp += (ll)num[a[l]] * num[a[l]];
70             l++;
71         }
72         while (l > interval[i].l)
73         {
74             l--;
75             tmp -= (ll)num[a[l]] * num[a[l]];
76             num[a[l]]++;
77             tmp += (ll)num[a[l]] * num[a[l]];
78         }
79         ans[interval[i].id].a = tmp - (r - l + 1);
80         ans[interval[i].id].b = (ll)(r - l + 1) * (r - l);
81         ans[interval[i].id].reduce();
82     }
83 }
84 
85 int main()
86 {
87     scanf("%d%d", &n, &m);
88     rep1(i, 1, n) scanf("%d", &a[i]);
89     rep1(i, 1, m)
90     {
91         interval[i].id = i; scanf("%d%d", &interval[i].l, &interval[i].r);
92     }
93     unit = (int)sqrt(n);
94     sot(interval, m);
95     solve();
96     rep1(i, 1, m) printf("%lld/%lld\n", ans[i].a, ans[i].b);
97     return 0;
98 }

View Code

曼哈顿MST写法：

  1 #include <cstdio>
  2 #include <cstdlib>
  3 #include <algorithm>
  4 #define N 50000
  5 #define Q 50000
  6 #define INFI 123456789
  7 
  8 typedef long long ll;
  9 struct edge
 10 {
 11     int next, node;
 12 } e[Q << 1 | 1];
 13 int head[N + 1], tot = 0;
 14 struct point
 15 {
 16     int x, y, n;
 17     bool operator < (const point &p) const
 18     {
 19         return x == p.x ? y < p.y : x < p.x;
 20     }
 21 } interval[Q + 1], p[Q + 1];
 22 struct inedge
 23 {
 24     int a, b, w;
 25     bool operator < (const inedge &x) const
 26     {
 27         return w < x.w;
 28     }
 29 } ie[Q << 3 | 1];
 30 int cnt = 0;
 31 struct BITnode
 32 {
 33     int w, p;
 34 } arr[Q + 1];
 35 int n, q, col[N + 1], a[Q + 1], *l[Q + 1], f[N + 1], c[N + 1];
 36 ll cur, ans[Q + 1];
 37 bool v[Q + 1];
 38 
 39 template <typename T>
 40 inline T abs(T x)
 41 {
 42     return x < (T)0 ? -x : x;
 43 }
 44 
 45 inline int dist(const point &a, const point &b)
 46 {
 47     return abs(a.x - b.x) + abs(a.y - b.y);
 48 }
 49 
 50 inline void addinedge(int a, int b, int w)
 51 {
 52     ++cnt;
 53     ie[cnt].a = a, ie[cnt].b = b, ie[cnt].w = w;
 54 }
 55 
 56 inline void addedge(int a, int b)
 57 {
 58     e[++tot].next = head[a];
 59     head[a] = tot, e[tot].node = b;
 60 }
 61 
 62 inline bool cmp(int *a, int *b)
 63 {
 64     return *a < *b;
 65 }
 66 
 67 inline int query(int x)
 68 {
 69     int r = INFI, p = -1;
 70     for (; x <= q; x += x & -x)
 71         if (arr[x].w < r) r = arr[x].w, p = arr[x].p;
 72     return p;
 73 }
 74 
 75 inline void modify(int x, int w, int p)
 76 {
 77     for (; x > 0; x -= x & -x)
 78         if (arr[x].w > w) arr[x].w = w, arr[x].p = p;
 79 }
 80 
 81 int find(int x)
 82 {
 83     return x == f[x] ? x : f[x] = find(f[x]);
 84 }
 85 
 86 inline ll calc(int x)
 87 {
 88     return (ll)x * (x - 1);
 89 }
 90 
 91 inline void add(int l, int r)
 92 {
 93     for (int i = l; i <= r; ++i)
 94     {
 95         cur -= calc(c[col[i]]);
 96         cur += calc(++c[col[i]]);
 97     }
 98 }
 99 
100 inline void remove(int l, int r)
101 {
102     for (int i = l; i <= r; ++i)
103     {
104         cur -= calc(c[col[i]]);
105         cur += calc(--c[col[i]]);
106     }
107 }
108 
109 void dfs(int x, int l, int r)
110 {
111     v[x] = true;
112     //Process right bound
113     if (r < interval[x].y)
114         add(r + 1, interval[x].y);
115     else if (r > interval[x].y)
116         remove(interval[x].y + 1, r);
117     //Process left bound
118     if (l < interval[x].x)
119         remove(l, interval[x].x - 1);
120     else if (l > interval[x].x)
121         add(interval[x].x, l - 1);
122     ans[x] = cur;
123     //Moving on to next query
124     for (int i = head[x]; i; i = e[i].next)
125     {
126         if (v[e[i].node]) continue;
127         dfs(e[i].node, interval[x].x, interval[x].y);
128     }
129     //Revert changes
130     //Process right bound
131     if (r < interval[x].y)
132         remove(r + 1, interval[x].y);
133     else if (r > interval[x].y)
134         add(interval[x].y + 1, r);
135     //Process left bound
136     if (l < interval[x].x)
137         add(l, interval[x].x - 1);
138     else if (l > interval[x].x)
139         remove(interval[x].x, l - 1);
140 }
141 
142 int main()
143 {
144     //Initialize
145     scanf("%d%d", &n, &q);
146     for (int i = 1; i <= n; ++i) scanf("%d", col + i);
147     for (int i = 1; i <= q; ++i) scanf("%d%d", &interval[i].x, &interval[i].y);
148     //Manhattan MST
149     for (int i = 1; i <= q; ++i) p[i] = interval[i], p[i].n = i;
150     for (int dir = 1; dir <= 4; ++dir)
151     {
152         //Coordinate transform
153         if (dir == 2 || dir == 4)
154             for (int i = 1; i <= q; ++i) std::swap(p[i].x, p[i].y);
155         else if (dir == 3)
156             for (int i = 1; i <= q; ++i) p[i].x = -p[i].x;
157         //Sort by x-coordinate
158         std::sort(p + 1, p + q + 1);
159         //Discretize
160         for (int i = 1; i <= q; ++i) a[i] = p[i].y - p[i].x, l[i] = &a[i];
161         std::sort(l + 1, l + q + 1, cmp);
162         int cnt = 1;
163         for (int i = 2; i <= q; ++i)
164             if (*l[i] == *l[i - 1]) *l[i - 1] = cnt;
165             else *l[i - 1] = cnt++;
166         *l[q] = cnt;
167         //Initialize BIT
168         for (int i = 1; i <= q; ++i) arr[i].w = INFI, arr[i].p = -1;
169         //Find points and add edges
170         for (int i = q; i > 0; --i)
171         {
172             int pos = query(a[i]);
173             if (pos != -1)
174                 addinedge(p[i].n, p[pos].n, dist(p[i], p[pos]));
175             modify(a[i], p[i].x + p[i].y, i);
176         }
177     }
178     //Kruskal
179     std::sort(ie + 1, ie + cnt + 1);
180     //Initialize disjoint set
181     for (int i = 1; i <= q; ++i) f[i] = i;
182     //Add edges
183     for (int i = 1; i <= cnt; ++i)
184         if (find(ie[i].a) != find(ie[i].b))
185         {
186             f[find(ie[i].a)] = find(ie[i].b);
187             addedge(ie[i].a, ie[i].b), addedge(ie[i].b, ie[i].a);
188         }
189 
190     //Modui Algorithm
191     ++c[col[1]];
192     dfs(1, 1, 1);
193     //Output
194     for (int i = 1; i <= q; ++i)
195     {
196         ll x = ans[i], y = calc(interval[i].y - interval[i].x + 1);
197         if (!x) printf("0/1\n");
198         else
199         {
200             ll g = __gcd(x, y);
201             printf("%lld/%lld\n", x / g, y / g);
202         }
203     }
204     return 0;
205 }

View Code

一些相关的题目：

bzoj 3289, 3236, 3585, 2120, 1878

SPOJ D-query

hdu 6278

zoj 4008

poj 2104 (然而这题肯定主席树更清晰)

Reference:

莫队算法 (Mo's Algorithm): https://zhuanlan.zhihu.com/p/25017840

莫队算法讲解（详尽版）: https://blog.csdn.net/thinfatty/article/details/72581276

曼哈顿距离最小生成树与莫队算法： https://blog.csdn.net/huzecong/article/details/8576908

posted @ 2019-05-15 21:05 JHSeng 阅读(1982) 评论(0) 编辑收藏举报

刷新页面返回顶部

JHSeng's Blog

必求甚解, 知难而进, 精益求精

莫队算法 (Mo's Algorithm)

公告