莫隊算法 (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 }