HDU-4609 3-idiots (FFT/NTT)

HDU-4609(FFT/NTT)

题意: 给出n个木棒，现从中不重复地选出3根来，求能拼出三角形的概率。

计算合法概率容易出现重复，所以建议计算不合法方案数

枚举选出的最大边是哪条，然后考虑剩下两条边之和小于等于它

两条边之和为\(x\)的方案数可以\(FFT/NTT\)得到，是一个简单的构造

即\(f(x)=\sum x^{length_i}\)，求出\(f(x)^2\)，就能得到和的方案数，但是会重复，包括自己和自己算，一对算两次

处理一下前缀和即可

#include<bits/stdc++.h>
using namespace std;

#define reg register
typedef long long ll;
#define rep(i,a,b) for(reg int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(reg int i=a,i##end=b;i>=i##end;--i)

template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); } 
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); } 

char IO;
int rd(){
	int s=0,f=0;
	while(!isdigit(IO=getchar())) if(IO=='-') f=1;
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

const double PI=acos(-1);
const int N=262150,P=998244353;
const int g=3;

int n,m;
int b[100010];
int rev[N];
typedef complex <double> Cp;
Cp d[N];
ll a[N];

void FFT(int n,Cp *a,int f){
	rep(i,0,n-1) if(rev[i]>i) swap(a[i],a[rev[i]]);
	for(reg int i=1;i<n;i<<=1) {
		Cp w(cos(PI/i),f*sin(PI/i));
		for(reg int l=0;l<n;l+=i*2) {
			Cp e(1,0);
			for(reg int j=l;j<l+i;j++,e=e*w) {
				Cp t=a[j+i]*e;
				a[j+i]=a[j]-t;
				a[j]=a[j]+t;
			}
		}
	}
}

int main(){
	rep(kase,1,rd()) {
		n=rd();
		int ma=0;
		rep(i,1,n) b[i]=rd(),d[b[i]]+=Cp(1,0),ma=max(ma,b[i]);
		int R=1,c=0;
		while(R<=ma*2) R<<=1,c++;
		rep(i,1,R) rev[i]=(rev[i>>1]>>1)|((i&1)<<(c-1));
		FFT(R,d,1);
		rep(i,0,R-1) d[i]=d[i]*d[i];
		FFT(R,d,-1);
		rep(i,0,R-1) a[i]=(ll)(d[i].real()/R+0.5);
		rep(i,1,n) a[b[i]*2]--;
		rep(i,1,R) a[i]/=2;    //处理重复情况
		rep(i,1,R) a[i]+=a[i-1];
		ll ans=0,sum=0;
		rep(i,1,n) {
			ans+=a[b[i]]; // 计算不合法方案数
			sum+=1ll*(i-1)*(i-2)/2;
		}
		rep(i,0,R+3) a[i]=0,d[i]=Cp(0,0);
		rep(i,0,ma+3) b[i]=0;
		ans=sum-ans;
		printf("%.7lf\n",1.0*ans/sum);
	}
}