【BZOJ 3513】idiots

传送门

---

Problem

给定 $n$ 个长度分别为 $a_i$ 的木棒，问随机选择 $3$ 个木棒能够拼成三角形的概率。

数据范围：$3\le n\le 10^5$，$1\le a_i\le 10^5$。

---

Solution

我们设 $cn{t}_i$ 表示 $i$ 的出现次数，用 $f_i$ 表示用两根木棒能拼出长度为 $i$ 的方案数，有：

$$f_i=\sum _{j=0}^{i}cn{t}_i\times cn{t}_{i-j}$$

这是一个卷积的式子，可以 $fft$ 优化。

要注意两点：

1. 当 $i$ 是偶数时，$f_i$ 就会加上 $cn{t}_{i/2}\times cn{t}_{i/2}$，就有可能会有选两根相同的木棒，这是不合法的。
2. 会算重，例如 $f_4$ 会算到两个 $cn{t}_1\times cn{t}_3$，因此最后要除 $2$。

因此对于 $i$，若 $i$ 是奇数， $F_i=\frac{f_i}{2}$；否则 $F_i=\frac{f_i-cn{t}_{i/2}}{2}$。

我们记后缀和 $S_i$ 表示长度大于等于 $i$ 的木棒有多少个。

那么不能拼出来的数量就是 $\sum F_i\times S_i$，由于总数是 $C_n^3=\frac{n!}{3!\times (n-3)!}=\frac{n(n-1)(n-2)}{6}$，答案就很显然了。

---

Code

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define N 300005
#define ll long long
using namespace std;
int n,Max,limit,pos[N];
const double pi=acos(-1);
struct complex{
	double x,y;
	complex(){}
	complex(double x,double y):x(x),y(y){}
	friend complex operator+(const complex &a,const complex &b)  {return complex(a.x+b.x,a.y+b.y);}
	friend complex operator-(const complex &a,const complex &b)  {return complex(a.x-b.x,a.y-b.y);}
	friend complex operator*(const complex &a,const complex &b)  {return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
	friend complex operator/(const complex &a,double b)  {return complex(a.x/b,a.y/b);}
}a[N];
ll A,B,f[N];
int g[N],cnt[N];
void Clear(){
	Max=0,limit=1;
	memset(a,0,sizeof(a));
	memset(g,0,sizeof(g));
	memset(cnt,0,sizeof(cnt));
}
void DFT(complex f[],int type){
	for(int i=0;i<limit;++i)
		if(pos[i]>i)  swap(f[i],f[pos[i]]);
	for(int mid=1;mid<limit;mid<<=1){
		complex now=complex(cos(pi/mid),type*sin(pi/mid));
		for(int i=0;i<limit;i+=(mid<<1)){
			complex w=complex(1,0);
			for(int j=0;j<mid;++j,w=w*now){
				complex p0=f[i+j],p1=w*f[i+j+mid];
				f[i+j]=p0+p1,f[i+j+mid]=p0-p1;
			}
		}
	}
	if(type==-1)  for(int i=0;i<limit;++i)  f[i].x/=limit;
}
int main(){
	int T;
	scanf("%d",&T);
	while(T--){
		Clear();
		scanf("%d",&n);
		for(int i=1,x;i<=n;++i){
			scanf("%d",&x);
			a[x].x++,cnt[x]++;
			Max=max(Max,x);
		}
		for(int i=Max;i>=1;--i)  g[i]=g[i+1]+cnt[i];
		while(limit<=(Max<<1))  limit<<=1;
		for(int i=0;i<limit;++i)  pos[i]=(pos[i>>1]>>1)|((i&1)*(limit>>1));
		DFT(a,1);
		for(int i=0;i<limit;++i)  a[i]=a[i]*a[i];
		DFT(a,-1);
		for(int i=1;i<=Max;++i){
			f[i]=(ll)(a[i].x+0.5);
			if(!(i&1))  f[i]-=cnt[i>>1];
			f[i]>>=1;
		}
		A=B=(ll)n*(n-1)*(n-2)/6;
		for(int i=1;i<=Max;++i)  A-=(ll)f[i]*g[i];
		printf("%.7lf\n",(double)A/(double)B);
	}
	return 0;
}