莫队算法初步

莫队算法，搞定！

思路很简单：对于每个询问，我们把它看成(l[i],r[i])这样一个平面上的点，那么询问之间的转移可以看成点之间的曼哈顿距离。那么我们要找到最小的距离怎么做呢？
1.曼哈顿路径，这个不好搞。
2.曼哈顿最小生成树，据说这个的权和是1的2倍，总之可以搞。

那么我们先要搞定manhattan最小生成树。如有n个点，那么复杂度是n^2的。这样肯定不行。
根据定理，每个点出发向8个方向的范围内最多只与一个点相连，那么我们只要想这8个点连边建图，这样只会有8n条边（去重了还剩4n），把这个图建出来，最小生成树就很好搞了。

给个链接：

http://wenku.baidu.com/view/1e4878196bd97f192279e941.html

其实莫队的实质就是找到一条转移代价最小的计算询问的顺序，那么我们只要按着这个顺序dfs地计算询问即可。注：计算的时候，用根去更新儿子，进行增减区间，然后注意回溯。Pascal果断还是慢啊。

const
    maxn=500000;
type
    et=record
        t,next:longint;
    end;
    tt=record
        n,s:longint;
    end;
    arr=array[0..maxn]of longint;
var
    t:array[0..maxn]of tt;
    e:array[0..maxn]of et;
    f,fir,u,v,val,l,r,s,a,b,key,pos:arr;
    col,ans,size:array[0..maxn]of int64;
    vis:array[0..maxn]of boolean;
    i,n,m,tot,d,rot:Longint;
    num,sum:int64;
procedure swap(var a,b:longint);
var
    tmp:longint;
begin
    tmp:=a; a:=b; b:=tmp;
end;
procedure sort(var a,b,c:arr; l,r:longint);
var
    i,j,x,z:longint;
begin
    i:=l; j:=r;
    x:=a[(l+r)>>1];
    z:=b[(l+r)>>1];
    repeat
        while (a[i]<x)or((a[i]=x)and(b[i]<z)) do inc(i);
        while (x<a[j])or((x=a[j])and(z<b[j])) do dec(j);
        if not (i>j) then
            begin
                swap(a[i],a[j]);
                swap(b[i],b[j]);
                swap(c[i],c[j]);
                inc(i); dec(j);
            end;
    until i>j;
    if l<j then sort(a,b,c,l,j);
    if i<r then sort(a,b,c,i,r);
end;

function find(i:longint):longint;
begin
    if f[i]=0 then exit(i);
    if f[f[i]]=0 then exit(f[i]);
    find:=find(f[i]);
    f[i]:=find;
end;

procedure add(i,x,k:longint);
begin
    while i>0 do
        begin
            if x<t[i].s then
                begin
                    t[i].s:=x; t[i].n:=k;
                end;
            i:=i-(i and -i);
        end;
end;
function query(i:longint):longint;
var
    tmp:longint;
begin
    query:=-1; tmp:=maxlongint;
    while i<=m do
        begin
            if t[i].s<tmp then
                begin
                    tmp:=t[i].s; query:=t[i].n;
                end;
            i:=i+(i and -i);
        end;
end;

procedure addedge(x,y,z:longint);
begin
    inc(tot);
    u[tot]:=x; v[tot]:=y; val[tot]:=z;
end;
procedure connect(x,y:longint);
begin
    inc(d);
    e[d].t:=y; e[d].next:=fir[x]; fir[x]:=d;
end;

procedure manhattan(x,y:arr);
var
    i,k,tmp,fx,fy:longint;
begin
    for k:=1 to 4 do
        begin
            if (k=2)or(k=4) then for i:=1 to m do swap(x[i],y[i])
            else if k=3 then for i:=1 to m do x[i]:=-x[i];
            sort(x,y,s,1,m);
            for i:=1 to m do key[i]:=y[i]-x[i];
            for i:=1 to m do pos[i]:=i;
            sort(key,pos,b,1,m);
            for i:=1 to m do key[pos[i]]:=i;
            for i:=1 to m do begin t[i].n:=-1; t[i].s:=maxlongint; end;
            for i:=m downto 1 do
                begin
                    tmp:=query(key[i]);
                    if tmp<>-1 then
                        addedge( s[i] , s[tmp] , abs(x[i]-x[tmp])+abs(y[i]-y[tmp]) );
                    add(key[i],x[i]+y[i],i);//这里我们插入的位置是key，比较的是x[i]+y[i]，返回的是i；
                end;
        end;
    sort(val,u,v,1,tot);
    for i:=1 to tot do
        begin
            fx:=find(u[i]); fy:=find(v[i]);
            if fx<>fy then
                begin
                    connect(u[i],v[i]);
                    connect(v[i],u[i]);
                    f[fx]:=v[i];
                end;
        end;
end;

procedure change(x,y:longint);
var
    i:longint;
begin
    ans[y]:=ans[x];
    for i:=l[x] to l[y]-1 do
        begin
            ans[y]:=ans[y]-2*(col[a[i]]-1);
            dec(col[a[i]]);
        end;
    for i:=r[y]+1 to r[x] do
        begin
            ans[y]:=ans[y]-2*(col[a[i]]-1);
            dec(col[a[i]]);
        end;
    for i:=l[y] to l[x]-1 do
        begin
            ans[y]:=ans[y]+2*col[a[i]];
            inc(col[a[i]]);
        end;
    for i:=r[x]+1 to r[y] do
        begin
            ans[y]:=ans[y]+2*col[a[i]];
            inc(col[a[i]]);
        end;
end;

procedure dfs(i:longint);
var
    j,k:longint;
begin
    vis[i]:=true;
    j:=fir[i];
    while j>0 do
        begin
            k:=e[j].t;
            if not vis[k] then
                begin
                    change(i,k);
                    dfs(k);
                    change(k,i);
                end;
            j:=e[j].next;
        end;
end;

function gcd(x,y:int64):int64;
begin
    if x=0 then exit(y)
    else exit(gcd(y mod x,x));
end;

begin
    //assign(input,'socks.in'); reset(input);
    //assign(output,'socks.out'); rewrite(output);
    readln(n,m);
    for i:=1 to n do read(a[i]);
    for i:=1 to m do
        begin
            readln(l[i],r[i]);
            if l[i]>r[i] then swap(l[i],r[i]);
            s[i]:=i;
            size[i]:=r[i]-l[i]+1;
        end;
    manhattan(l,r);
    randomize;
    rot:=random(m)+1;
    for i:=l[rot] to r[rot] do
        inc(col[a[i]]);
    for i:=1 to n do if col[i]>1 then ans[rot]:=ans[rot]+col[i]*(col[i]-1);
    dfs(rot);
    for i:=1 to m do
        begin
            if size[i]=1 then sum:=1
            else sum:=size[i]*(size[i]-1);
            num:=gcd(ans[i],sum);
            writeln(ans[i] div num ,'/',sum div num);
        end;
end.

 

---

 

Update:

　　我这个弱B，一直不知道莫队有一种更为简洁的写法，导致这篇博文在这挂了这么久，被坑到的同学是在对不起......

　　真是的，要是莫队算法这么长的话还有谁愿去写！！！

　　其实不用manhattan最小生成树了啦，直接把x轴分块，每一块中按y升序排列就好了。这样每个块内做一遍暴力的莫队，时间算下来O(n^1.5)。具体实现是按照x/s,y排序(x是根号m)，然后对于x/s相同的先把第一个求出来，然后暴力转移后面的。

　　这才是正宗的得分利器莫对算法嘛！

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#define maxn 55000
#define inf 2147483647
using namespace std;
struct query
{
    int l,r,s,w;
}a[maxn];
int c[maxn];
long long col[maxn],size[maxn],ans[maxn];
int n,m,cnt,len;

long long gcd(long long x,long long  y)
{
    return (!x)?y:gcd(y%x,x);
}

bool cmp(query a,query b)
{
    return (a.w==b.w)?a.r<b.r:a.w<b.w;
}

int main()
{
    //freopen("hose.in","r",stdin);
    scanf("%d%d",&n,&m);
    for (int i=1;i<=n;i++) scanf("%d",&c[i]);
    len=(int)sqrt(m);
    cnt=(len*len==m)?len:len+1;
    for (int i=1;i<=m;i++) 
    {
        scanf("%d%d",&a[i].l,&a[i].r);
        if (a[i].l>a[i].r) swap(a[i].l,a[i].r);
        size[i]=a[i].r-a[i].l+1;
        a[i].w=a[i].l/len+1;
        a[i].s=i;
    }
    sort(a+1,a+m+1,cmp);
    int i=1;
    while (i<=m)
    {
        int now=a[i].w;
        memset(col,0,sizeof(col));
        for (int j=a[i].l;j<=a[i].r;j++) ans[a[i].s]+=2*(col[c[j]]++);
        i++;
        for (;a[i].w==now;i++)
        {
            ans[a[i].s]=ans[a[i-1].s];
            for (int j=a[i-1].r+1;j<=a[i].r;j++)
                ans[a[i].s]+=2*(col[c[j]]++);
            if (a[i-1].l<a[i].l)
                for (int j=a[i-1].l;j<a[i].l;j++) 
                    ans[a[i].s]-=2*(--col[c[j]]);
            else
                for (int j=a[i].l;j<a[i-1].l;j++)
                    ans[a[i].s]+=2*(col[c[j]]++);
        }
    }
    long long all,num;
    for (int i=1;i<=m;i++)
    {
        if (size[i]==1) all=1; else all=size[i]*(size[i]-1);
        num=gcd(ans[i],all);
        printf("%lld/%lld\n",ans[i]/num,all/num);
    }
    return 0;
}