莫队算法及其应用

在写这篇博客之前，我最想做的一件事就是：ORZ莫队%%%%%%%%。

说明：ceil(x)表示x向上取整，sqrt(x)表示对x开算数平方根。

一、莫队算法简介

　　莫队算法是一种暴力算法，真的很暴力，但速度很快，属于速度快的暴力。它的基本思想就是分块。关于分块的介绍建议参考hzwer的博客，然后%%%%hzw。莫队算法主要用于解决一类离线查询的问题，和线段树处理的问题是一样的，但处理的是两个不同的方面，当由[L,R]转移到[L’，R’]的时间为O(|L‘-L|+|R‘-R|)时适宜使用莫队算法。这个可以从题目中体会。因为采取的是分块，它的复杂度是O(nsqrt(n))。其实质是将询问按照某种顺序排好，这个也应该从题目中去体会，我们参考一道题目。

二、典型例题

　　著名例题，小Z的袜子

　　链接：http://www.lydsy.com/JudgeOnline/problem.php?id=2038

　　题目是中文的，看得懂所以不复制粘贴了。题意也不难理解，稍有组合数学常识的人都可以看出。

三、解法

　　因为题目中的组合是C(n,2)，所以我们预处理出C2数组，存放2-n对2的组合数，作为特例，C(0，2)=C(1，2)=0;

　　我们用桶tab存放[L,R]中每种颜色的数量，假设我们求出了[L,R]，求[L+1,R](或[L-1,R][L,R+1][L,r-1])时只需要把桶里面的--或++就可以了，令[L,R]的答案为ans，那么[L+1,R]的答案为ans-C(tab[L],2)+C(tab[L]+1,2),这是O(1)的;

　　我们可以发现，假设我们求出了[L,R]，那么我们求出[L’，R’]的时间为O(|L‘-L|+|R‘-R|)，所以我们采用莫队算法。

　　数据范围是n，m<=50000，这启发我们用分块（当然如果执意要写曼哈顿最小生成树那也没人拦你）。我们先将所有询问按照l为第一关键字，r为第二关键字排一遍序，再将排好序的数组分成[√n]块,再将分好块的数组按照r大小排一遍序，这样我们就做完了第一步了。

　　接着我们按块处理，对于每一块，找出每个询问和它前面一个询问的差异，修改差异，不断地这么做，就可以得到答案。

　　这样做总时间复杂度仅有O(n√n),比原有的O(n^2)的暴力快了许多，但这是为什么呢？

四、复杂度分析

　　首先是分块这一步，这一步的时间复杂度毫无疑问地是O(√n*√n*log√n+nlogn)=O(nlogn);

　　接着就到了莫队算法的精髓了，下面我们用通俗易懂的初中方法来证明它的时间复杂度是O(n√n)；

　　证：令每一块中L的最大值为max₁,max₂,max₃,...,max_ceil(√n).

　　由第一次排序可知，max₁<=max₂<=...<=max_ceil(√n)

　　显然，对于每一块暴力求出第一个询问的时间复杂度为O(n)。

　　考虑最坏的情况，在每一块中，R的最大值均为n，每次修改操作均要将L由max_i-1修改至max_i或由max_i修改至max_i-1。

　　考虑R：因为R在块中已经排好序，所以在同一块修改完它的时间复杂度为O(n)。对于所有块就是O(n√n)。

　　重点分析L：因为每一次改变的时间复杂度都是O(max_i-max_i-1)的，所以在同一块中时间复杂度为O(√n*(max_i-max_i-1)).

　   将每一块L的时间复杂度合在一起，可以得到对于L的总时间复杂度为

　　　　O(√n*(max₁-1)+√n*(max₂-max₁)+√n*(max₃-max₂)+...+√n*(max_ceil(√n)-max_ceil(√n-1)))

　　    =O(√n*(max₁-1+max₂-max₁+max₃-max₂+...+max_ceil(√n-1)-max_ceil(√n-2)+max_ceil(√n)-max_ceil(√n-1)))

　　    =O(√n*(max_ceil(√n)-1))  (初中裂项求和)

　　由题可知max_ceil(√n)最大为n，所以L的总时间复杂度最坏情况下为O(n√n).

　　综上所述，莫队算法的时间复杂度为O(n√n)；

五、例题代码

　　还是用emacs写的，所以还是两格缩进，不喜勿喷。

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 typedef long long ll;
 4 ll a[60000],tab[60000];
 5 struct ask{
 6   ll l,r,num;
 7 }b[60000];
 8 ll cmp(ask x,ask y){
 9   if(x.l<y.l) return 1;
10   if(x.l>y.l) return 0;
11   if(x.r<y.r) return 1;
12   return 0;
13 }
14 ll comp(ask x,ask y){
15   if(x.r<y.r) return 1;
16   if(x.r>y.r) return 0;
17   if(x.l<y.l) return 1;
18   return 0;
19 }
20 ll gcd(ll a,ll b){
21   if(!b) return a;
22   return gcd(b,a%b);
23 }ll n,m;
24 ll comb2[60000];//组合数C(n,2)
25 ll prix[60000],priy[60000];//答案
26 ll rep(ll ol,ll nl,ll lr,ll nr,ll &ans){//回答修改的问题，原来的是[ol,lr],现在是[nl,nr];
27   if(ol<=nl)
28     for(ll i=ol;i<nl;i++){ans-=comb2[tab[a[i]]]; tab[a[i]]--;ans+=comb2[tab[a[i]]];}
29   else
30     for(ll i=ol-1;i>=nl;i--){ans-=comb2[tab[a[i]]]; tab[a[i]]++;ans+=comb2[tab[a[i]]];}
31   for(ll i=lr+1;i<=nr;i++){ans-=comb2[tab[a[i]]]; tab[a[i]]++;ans+=comb2[tab[a[i]]];}
32   return ans;
33 }
34
35 int main(){
36   scanf("%lld%lld",&n,&m);comb2[1]=comb2[0]=0;
37   for(ll i=2;i<=n;i++)comb2[i]=(ll)((double)i/2.0*(double)(i-1));//计算组合数
38   for(ll i=1;i<=n;i++) scanf("%lld",&a[i]);
39   for(ll i=1;i<=m;i++){
40     scanf("%lld%lld",&b[i].l,&b[i].r);
41     b[i].num=i;
42   }
43   ll sq=sqrt(m);
44   sort(b+1,b+m+1,cmp);//第一次排序
45   for(ll i=1;i<=m;i+=sq){
46     sort(b+i,b+min(i+sq,m+1),comp);//第二次排序
47   }
48   for(ll i=1;i<=m;i+=sq){
49     ll ed=min(m,i+sq-1);
50     memset(tab,0,sizeof(tab));ll maxx=0;
51     long long ans=0;ans=rep(b[i].l,b[i].l,b[i].l-1,b[i].r,ans);//同下
52     prix[b[i].num]=ans;priy[b[i].num]=comb2[b[i].r-b[i].l+1];//暴力算出每块的第一个，其实这里可以不这么做，直接继承上一块也行
53     if(prix[b[i].num]==0)priy[b[i].num]=1;
54     else{ll g=gcd(prix[b[i].num],priy[b[i].num]);
55       prix[b[i].num]/=g;priy[b[i].num]/=g;}//约分
56     for(ll j=i+1;j<=ed;j++){
57       prix[b[j].num]=rep(b[j-1].l,b[j].l,b[j-1].r,b[j].r,ans);//从上一个询问推导这一个询问
58       priy[b[j].num]=comb2[b[j].r-b[j].l+1];
59       if(prix[b[j].num]==0)priy[b[j].num]=1;
60       else{
61    　　 ll g=gcd(prix[b[j].num],priy[b[j].num]);
62     　　prix[b[j].num]/=g;priy[b[j].num]/=g;
63       }
64     }
65   }
66   for(ll i=1;i<=m;i++){
67     printf("%lld/%lld\n",prix[i],priy[i]);//这里需要注意，BZOJ有坑，cout是会RE的
68   }
69   return 0;
70 }