bzoj3513 [MUTC2013]idiots

Description

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

Input

第一行T(T<=20)，表示数据组数。 接下来若干行描述T组数据，每组数据第一行是n，接下来一行有n个数表示a_i。1<=a_i<=100000，N<=100000

Output

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

设a(x)表示长度x的木棒出现次数

b(x)表示选两根木棒（可重）长度之和为x的方案数

显然b(m)=Σ a(x)*a(m-x) x=0..m，可以用fft求

稍加处理即可得到合法方案数

#include<bits/stdc++.h>
int N,P;
const double pi=3.14159265358979323846;
int r[262144];
struct cplx{
    double a,b;
    inline cplx(double r=0,double i=0){a=r;b=i;}
    inline cplx operator+(cplx x){return cplx(a+x.a,b+x.b);}
    inline cplx operator-(cplx x){return cplx(a-x.a,b-x.b);}
    inline cplx operator*(cplx x){return cplx(a*x.a-b*x.b,a*x.b+b*x.a);}
}a[262144];
void dft(cplx*a,int t){
    for(int i=0;i<N;i++)if(i>r[i])std::swap(a[i],a[r[i]]);
    for(int i=1;i<N;i<<=1){
        cplx w(cos(pi/i),t*sin(pi/i));
        for(int j=0;j<N;j+=i<<1){
            cplx e(1,0),*b=a+j,*c=b+i;
            for(int k=0;k<i;k++,e=e*w){
                cplx x=b[k],y=c[k]*e;
                b[k]=x+y,c[k]=x-y;
            }
        }
    }
}
inline int input(){
    int x=0,c=getchar();
    while(c>57||c<48)c=getchar();
    while(c>47&&c<58)x=x*10+c-48,c=getchar();
    return x;
}
int t[262144];
double t2[262144];
int main(){
    int T=input();
    while(T--){
        int n=input();
        for(int i=0;i<=100000;i++)t[i]=0;
        int mx=2;
        for(int i=0;i<n;i++){
            int x=input();
            ++t[x];
            if(x>mx)mx=x;
        }
        for(N=2,P=0;N<=mx;N<<=1,++P);
        if(N<262144)N<<=1,++P;
        for(int i=1;i<N;i++)r[i]=r[i>>1]>>1|(i&1)<<P;
        for(int i=0;i<=mx;i++)a[i]=cplx(t[i],0);
        for(int i=mx+1;i<N;i++)a[i]=cplx(0,0);
        dft(a,1);
        for(int i=0;i<N;i++)a[i]=a[i]*a[i];
        dft(a,-1);
        for(int i=0;i<=mx;i++)t2[i]=a[i].a/N;
        for(int i=0;i<=mx>>1;i++)t2[i<<1]-=t[i];
        for(int i=1;i<=mx;i++)t2[i]+=t2[i-1];
        double c=0;
        for(int i=1;i<=mx;i++)c+=t2[i]*t[i];
        printf("%.7lf\n",1.0-c/(double(n)*(n-1)*(n-2)/3.0));
    }
    return 0;
}