莫队算法 (Mo's Algorithm)

时间 2020-06-20

标签算法 mo's algorithm 繁體版

原文原文链接

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

好比GDCPC2019广东省赛就有这么一道题：c++

　　给定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作法。数组

那应该怎么作呢？因为这道题没有修改操做，因此能够考虑离线处理，最明显的处理区间问题的离线算法就是莫队了……莫队本质上是一种分块算法，一般用来离线解决只查询不修改的区间问题。ide

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

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

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

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

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

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

分块相同时，右端点递增是的，分块共有 $O(\sqrt{N} )$ 个，复杂度为 $O(N^{1.5} )$
分块转移时，右端点最多变化，分块共有 $O(\sqrt{N} )$ 个，复杂度为 $O(N^{1.5} )$
分块相同时，左端点最多变化 $\sqrt{N}$ ，分块转移时，左端点最多变化 $2\sqrt{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

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

设在某一区间内共有颜色，每双袜子的个数为。

答案为 $(\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