HDU 5751 Eades（单调栈+FFT）

Description
Peter有一个序列a[1],a[2],…,a[n]，定义g(l,r)表示子序列{a[l],a[l+1],…,a[r]}的最大值, f(l,r)=sum{[a[i]==g(l,r)]}(l<=i<=r).
注意[condition]=1当且仅当condition是true, 否则[condition]=0
对于每个整数k∈{1,2,…,n}, Peter想要知道有多少整数对l和r (l≤r)满足f(l,r)=k
Input
输入包含多组数据, 第一行包含一个整数T表示测试数据组数. 对于每组数据:
第一行包含一个整数n(1≤n≤60000)表示序列的长度，第二行包含n个整数ai(1≤ai≤n)
Output
对于每组数据, 输出一个整数 ，其中z​[k]表示满足f(l,r)=k的数对l和r的个数, ⊕是异或位运算操作.
Sample Input
3
3
1 2 3
4
1 1 1 1
6
1 2 2 1 1 2
Sample Output
12
12
36
Solution
首先用单调栈预处理出每个数x作为最大值的最大区间[lx,rx]，假设在这个区间中x出现m次，出现位置分别为p[1],p[2],…,p[m]，那么这m个x就可以转化为一个长度为m+1的序列b，b[0]=p[1]-lx+1，b[m]=rx-p[m]+1，b[i]=p[i+1]-p[i]，i=1,2,…,m-1，简单推导可知这m个x对z[k]的贡献为，令x[i]=b[i],y[i]=b[m-i]，那么有，后者直接计算复杂度为O(m^2)，用fft对x[0],…,x[m]和y[0],…,y[m]做一遍卷积可以将复杂度降为O(mlogm)，注意此处计算的是m个x对答案的贡献
Code

#include<cstdio>#include<iostream>#include<cstring>#include<cmath>#include<vector>using namespace std;typedef long long ll;#define maxn 244444#define PI acos(-1.0) struct complex{    double r,i;    complex(double _r=0,double _i=0)    {        r=_r,i=_i;    }    complex operator +(const complex &b)    {        return complex(r+b.r,i+b.i);    }    complex operator -(const complex &b)    {        return complex(r-b.r,i-b.i);    }    complex operator *(const complex &b)    {        return complex(r*b.r-i*b.i,r*b.i+i*b.r);    }    complex conj()    {        return complex(r,-i);    }};int T,n,a[maxn],l[maxn],r[maxn],flag[maxn],b[maxn];vector<int>v[maxn];ll z[maxn],ans; complex x[maxn],y[maxn];int pos[maxn];void fft_init(int len){    int j=0;    while((1<<j)<len)j++;    j--;    for(int i=0;i<len;i++)        pos[i]=pos[i>>1]>>1|((i&1)<<j);}void fft(complex *x,int len,int sta){    for(int i=0;i<len;i++)        if(i<pos[i])swap(x[i],x[pos[i]]);    for(int m=2;m<=len;m<<=1)    {        complex Wn(cos(sta*2*PI/m),sin(sta*2*PI/m));        for(int i=0;i<len;i+=m)        {            complex W(1,0);            for(int j=i;j<i+m/2;j++)            {                complex x1=x[j],x2=W*x[j+m/2];                x[j]=x1+x2,x[j+m/2]=x1-x2;                W=W*Wn;            }        }    }    if(sta==-1)        for(int i=0;i<len;i++)            x[i].r/=len;}void init(){    ans=0;    memset(z,0,sizeof(z));    memset(flag,0,sizeof(flag));    for(int i=1;i<=n;i++)v[i].clear();}void deal(int *a,int n){    int p,sta[maxn];    p=0;    for(int i=1;i<=n;i++)    {        if(!p||a[i]<=a[sta[p]])sta[++p]=i;        else         {            while(p&&a[i]>a[sta[p]])                r[sta[p]]=i-1,p--;            sta[++p]=i;        }    }    while(p)r[sta[p]]=n,p--;    for(int i=n;i>=1;i--)    {        if(!p||a[i]<=a[sta[p]])sta[++p]=i;        else         {            while(p&&a[i]>a[sta[p]])                l[sta[p]]=i+1,p--;            sta[++p]=i;        }    }    while(p)l[sta[p]]=1,p--;}void cal(int *b,int m){    int len=1;    while(len<=2*m)len<<=1;    fft_init(len);    for(int i=0;i<=m;i++)x[i]=complex(b[i],b[m-i]);    for(int i=m+1;i<len;i++)x[i]=complex(0,0);    fft(x,len,1);    for(int i=0;i<len;i++)    {        int j=(len-i)&(len-1);        y[i]=(x[i]*x[i]-(x[j]*x[j]).conj())*complex(0,-0.25);    }    fft(y,len,-1);    for(int i=0;i<m;i++)    {        ll temp=(ll)(y[i].r+0.5);        z[m-i]+=temp;    }}int main(){    scanf("%d",&T);    while(T--)    {        init();        scanf("%d",&n);        for(int i=1;i<=n;i++)            scanf("%d",&a[i]),v[a[i]].push_back(i);        deal(a,n);        for(int i=1;i<=n;i++)            if(!flag[i])            {                int m=0,last=l[i]-1,j;                for(j=0;j<v[a[i]].size();j++)                {                    if(v[a[i]][j]>r[i])break;                    int now=v[a[i]][j];                    flag[now]=1;                    b[m++]=now-last;                    last=now;                }                while(v[a[i]].size()&&v[a[i]][0]<=r[i])                    v[a[i]].erase(v[a[i]].begin());                b[m++]=r[i]-last+1;                cal(b,m-1);            }        for(int i=1;i<=n;i++)ans+=i^z[i];        printf("%I64d\n",ans);    }    return 0;}