FFT+計數--bzoj3513: [MUTC2013]idiots
阿新 • • 發佈:2018-12-02
傳送門
因為不能構成三角形的三邊滿足
用
求出
的方案數,再把重複取一個的減掉,然後從小到大列舉第三條邊,把前面小於的都加上,再乘以當前長度的邊的個數加到答案裡
因為 和 是一樣的,所以要除以2
再用總數 減掉不合法的再除以總數就是概率
注意 ,寫成 貢獻數次
#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;
}