【BZOJ3509】【CodeChef】—COUNTARI（分块+FFT）

传送门

$BZOJ$看不了权限题怎么办？
把网址中的$problem$改成$show$就可以了(虽然没法提交)
可以从$darkbzoj$上下数据啊

考虑处理$2$个东西
$l[i][j]$表示$i$前面，值为$j$的个数，$r[i][j]$表示$i$后面，值为$j$的个数

则$ans=sum_{i=1}^{n}sum_{j=0}^{max(a)}l[i][j]*r[i][2*a[i]-j]$

发现这是一个卷积的形式，但直接$fft$也不行

考虑分块对于当前块统计这个块前的$l$，块后的$r$
对于有在块内的再重新统计

块大小大约在$sqrt {nlogn}$最佳
复杂度$O(nsqrt{nlogn})$

#include<bits/stdc++.h>
using namespace std;
#define ll long long
inline int read(){
	char ch=getchar();
	int res=0,f=1;
	while(!isdigit(ch)){if(ch=='-')f=-f;ch=getchar();}
	while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=getchar();
	return res*f;
}
const int N=100005;
const int M=120005;
const double pi=acos(-1);
struct plx{
	double x,y;
	plx(double _x=0,double _y=0):x(_x),y(_y){}
	friend inline plx operator +(const plx &a,const plx &b){
		return plx(a.x+b.x,a.y+b.y);
	}
	friend inline plx operator -(const plx &a,const plx &b){
		return plx(a.x-b.x,a.y-b.y);
	}
	friend inline plx operator *(const plx &a,const plx &b){
		return plx(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
	}
}A[N],B[N];
int rev[M],lim=1,tim=0;
inline void fft(plx f[],int kd){
	for(int i=0;i<lim;i++)if(i<rev[i])swap(f[i],f[rev[i]]);
	for(int mid=1;mid<lim;mid<<=1){
		plx now=plx(cos(pi/mid),kd*sin(pi/mid));
		for(int i=0;i<lim;i+=(mid<<1)){
			plx w=plx(1,0);
			for(int j=0;j<mid;j++,w=w*now){
				plx a0=f[i+j],a1=w*f[i+j+mid];
				f[i+j]=a0+a1,f[i+j+mid]=a0-a1;
			}
		}
	}
	if(kd==-1)for(int i=0;i<lim;i++)f[i].x/=lim;
}
int n,a[N],L[N],R[N],mx,cnt,blo,l[M],r[M];
ll ans;
int main(){
	n=read(),blo=2000,cnt=(n-1)/blo+1;
	for(int i=1;i<=n;i++)a[i]=read(),r[a[i]]++,mx=max(mx,a[i]);
	while(lim<=mx*2)lim<<=1,tim++;
	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(tim-1));
	for(int i=1;i<=cnt;i++)L[i]=(i-1)*blo+1,R[i]=i*blo;R[cnt]=n;
	for(int i=1;i<=cnt;i++){
		for(int j=L[i];j<=R[i];j++)r[a[j]]--;
		for(int j=0;j<lim;j++)A[j]=B[j]=plx(0,0);
		for(int j=1;j<=mx;j++)A[j]=plx(l[j],0),B[j]=plx(r[j],0);
		fft(A,1),fft(B,1);
		for(int j=0;j<lim;j++)A[j]=A[j]*B[j];
		fft(A,-1);
		for(int j=L[i];j<=R[i];j++)ans+=(ll)(A[2*a[j]].x+0.5);
		for(int j=L[i];j<=R[i];j++){
			for(int k=L[i];k<j;k++)
				if(2*a[j]-a[k]>0)ans+=r[2*a[j]-a[k]];
			for(int k=j+1;k<=R[i];k++)
				if(2*a[j]-a[k]>0)ans+=l[2*a[j]-a[k]];
			l[a[j]]++;
		}
	}
	cout<<ans<<'
';
}