【刷题】BZOJ 3513 [MUTC2013]idiots

Description

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

Input

第一行T(T<=100)，表示数据组数。

接下来若干行描述T组数据，每组数据第一行是n，接下来一行有n个数表示a_i。

3≤N≤10^{5,1≤a_i≤10}5

Output

T行，每行一个整数，四舍五入保留7位小数。

Sample Input

2
4
1 3 3 4
4
2 3 3 4

Sample Output

0.5000000
1.0000000

HINT

T<=20

N<=100000

Solution

首先求 \(b_i\) 表示一对木棒加起来的和不超过 \(i\) ，有多少种方案
这个用FFT先求正好是 \(i\) 长度有多少种方案，然后累前缀和就可以了
主要把一个木棒用两次的和只是换了个位置，实际是一样的方案去掉
接着，考虑每个木棒， \(b_{a_i}\) 的大小就是无法组成三角形的方案数
将它们求和
再用总方案数减去不合法方案数，就是合法方案数了
除一下就是答案

#include<bits/stdc++.h>
#define ui unsigned int
#define ll long long
#define db double
#define ld long double
#define ull unsigned long long
const int MAXN=1<<19;
const db Pi=acos(-1.0);
int T,ln,n,m,qn,rev[MAXN],a[MAXN],cnt;
ll sum,all,b[MAXN];
struct Complex{
	db real,imag;
	inline Complex operator + (const Complex &A) const {
		return (Complex){real+A.real,imag+A.imag};
	};
	inline Complex operator - (const Complex &A) const {
		return (Complex){real-A.real,imag-A.imag};
	};
	inline Complex operator * (const Complex &A) const {
		return (Complex){real*A.real-imag*A.imag,imag*A.real+real*A.imag};
	};
};
Complex x[MAXN];
template<typename T> inline void read(T &x)
{
	T data=0,w=1;
	char ch=0;
	while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar();
	if(ch=='-')w=-1,ch=getchar();
	while(ch>='0'&&ch<='9')data=((T)data<<3)+((T)data<<1)+(ch^'0'),ch=getchar();
	x=data*w;
}
template<typename T> inline void write(T x,char ch='\0')
{
	if(x<0)putchar('-'),x=-x;
	if(x>9)write(x/10);
	putchar(x%10+'0');
	if(ch!='\0')putchar(ch);
}
template<typename T> inline void chkmin(T &x,T y){x=(y<x?y:x);}
template<typename T> inline void chkmax(T &x,T y){x=(y>x?y:x);}
template<typename T> inline T min(T x,T y){return x<y?x:y;}
template<typename T> inline T max(T x,T y){return x>y?x:y;}
inline void FFT(Complex *A,int tp)
{
	for(register int i=0;i<n;++i)
		if(i<rev[i])std::swap(A[i],A[rev[i]]);
	for(register int l=2;l<=n;l<<=1)
	{
		Complex wn=(Complex){cos(2*Pi/l),sin(tp*2*Pi/l)};
		for(register int i=0;i<n;i+=l)
		{
			Complex w=(Complex){1,0};
			for(register int j=0;j<(l>>1);++j)
			{
				Complex A1=A[i+j],A2=w*A[i+j+(l>>1)];
				A[i+j]=A1+A2,A[i+j+(l>>1)]=A1-A2;
				w=w*wn;
			}
		}
	}
}
int main()
{
	read(T);
	while(T--)
	{
		read(ln);
		qn=0;
		for(register int i=1;i<=ln;++i)read(a[i]),chkmax(qn,a[i]);
		qn++;m=qn+qn-1;
		cnt=0;
		for(n=1;n<m;n<<=1)cnt++;
		for(register int i=0;i<n;++i)x[i].real=x[i].imag=0;
		for(register int i=1;i<=ln;++i)x[a[i]].real+=1.0;
		for(register int i=0;i<n;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(cnt-1));
		FFT(x,1);
		for(register int i=0;i<n;++i)x[i]=x[i]*x[i];
		FFT(x,-1);
		for(register int i=0;i<n;++i)b[i]=(int)(x[i].real/n+0.5);
		for(register int i=1;i<=ln;++i)b[a[i]<<1]--;
		for(register int i=0;i<n;++i)b[i]=(i?b[i-1]:0)+(b[i]>>1);
		all=1ll*ln*(ln-1)*(ln-2)/6;
		sum=0;
		for(register int i=1;i<=ln;++i)sum+=b[a[i]];
		printf("%.7f\n",(db)(all-sum)/all);
	}
	return 0;
}