FFT+计数--bzoj3513: [MUTC2013]idiots

传送门
因为不能构成三角形的三边满足$a+b\le c$
用$FFT$求出$a+b=c$的方案数，再把重复取一个的减掉，然后从小到大枚举第三条边，把前面小于的都加上，再乘以当前长度的边的个数加到答案里

因为$a,b$和$b,a$是一样的，所以要除以2

再用总数$C_n^3$减掉不合法的再除以总数就是概率

注意$limit<=2\ast n$，写成$n$贡献数次$wa$

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define maxn 400005
#define LL long long
using namespace std;
const double Pi=acos(-1.0);

inline int rd(){
	int x=0,f=1;char c=' ';
	while(c<'0' || c>'9') f=c=='-'?-1:1,c=getchar();
	while(c<='9' &&  c>='0') x=x*10+c-'0',c=getchar();
	return x*f;
}
inline int max(int x,int y){return x>y?x:y;}

int t,n,limit=1,l,rev[maxn],val[maxn],w[maxn],cnt[maxn];

struct complex{
	double x,y;
	complex(double xx=0,double yy=0){x=xx,y=yy;}
	void clear(){x=0,y=0;}
}a[maxn];
complex operator +(complex a,complex b){return complex(a.x+b.x,a.y+b.y);}
complex operator -(complex a,complex b){return complex(a.x-b.x,a.y-b.y);}
complex operator *(complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

inline void FFT(complex *F,int type){
	for(int i=0;i<limit;i++)
		if(i<rev[i]) swap(F[i],F[rev[i]]);
	for(int mid=1;mid<limit;mid<<=1){
		complex Wn(cos(Pi/mid),1.0*type*sin(Pi/mid));
		for(int r=mid<<1,j=0;j<limit;j+=r){
			complex w(1,0);
			for(int k=0;k<mid;k++,w=w*Wn){
				complex x=F[j+k],y=w*F[j+mid+k];
				F[j+k]=x+y,F[j+mid+k]=x-y;
			}
		}
	}
	if(type==-1) {
		for(int i=0;i<limit;i++) F[i].x/=limit;
	}
}

int main(){
	t=rd();
	while(t--){
		n=rd(); int mx=0; memset(cnt,0,sizeof cnt);
		for(int i=1;i<=n;i++) {
			w[i]=rd(); mx=max(mx,w[i]);
			cnt[w[i]]++;
		}
		limit=1; l=0;
		while(limit<=2*mx) limit<<=1,++l;
		for(int i=0;i<limit;i++)
			rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1)),a[i].clear();
		for(int i=1;i<=mx;i++) a[i].x=cnt[i];
		FFT(a,1);
		for(int i=0;i<limit;i++) a[i]=a[i]*a[i];
		FFT(a,-1);
		for(int i=1;i<=mx;i++) val[i]=(int)(a[i].x+0.5);
		for(int i=1;i<=mx;i++)
			if(!(i&1)) val[i]-=cnt[i/2];
		LL tmp=0,pre=0;
		for(int i=1;i<=mx;i++){
			pre+=val[i];
			tmp+=1LL*pre*cnt[i];
		}
		LL res=1LL*n*(n-1)*(n-2)/6;
		printf("%.7lf\n",1.0*(res-tmp/2)/res);
	}
	return 0;
}