bzoj 3513: [MUTC2013]idiots -- FFT

3513: [MUTC2013]idiots

Time Limit: 20 Sec  Memory Limit: 128 MB

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

Source

好像只能手写complex。。\好生气(︶︿︶)

不知道为什么我的常数巨大。。。

原谅我粘个题解qaq

#include<map>
#include<cmath>
#include<queue> 
#include<cstdio>
#include<complex>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
//#define cp complex<double>
#define inf 1000000007
#define ll long long
#define PI 3.1415926535897932384626433832795028841971693993751
#define N 400010
char xB[1<<15],*xS=xB,*xTT=xB;
#define getc() (xS==xTT&&(xTT=(xS=xB)+fread(xB,1,1<<15,stdin),xS==xTT)?0:*xS++)
#define isd(c) (c>='0'&&c<='9')
inline int rd(){
    char xchh;
    int xaa;
    while(xchh=getc(),!isd(xchh));(xaa=xchh-'0');
    while(xchh=getc(),isd(xchh))xaa=xaa*10+xchh-'0';return xaa;
}
struct cp{
     double r,i;
    cp(double _r=0,double _i=0):r(_r),i(_i){}
    cp operator + (cp x){return cp(r+x.r,i+x.i);}
    cp operator - (cp x){return cp(r-x.r,i-x.i);}
    cp operator * (cp x){return cp(r*x.r-i*x.i,r*x.i+i*x.r);}
    void clear(){r=i=0.0;}
};
cp a[N],b[N];
int r[N],m,L,M;
void fft(cp *x,int f)
{
    for(int i=0;i<m;i++) if(i<r[i]) swap(x[i],x[r[i]]);
    for(int i=1;i<m;i<<=1)
    {
        cp wn(cos(PI/i),f*sin(PI/i));
        for(int j=0;j<m;j+=(i<<1))
        {
            cp w(1,0),X,Y;
            for(int k=0;k<i;k++,w=w*wn)
            {
                X=x[j+k];Y=w*x[k+j+i];
                x[j+k]=X+Y;x[j+k+i]=X-Y;
            }
        }
    }
}
int T,n,c,p[N],f[N],mx;
ll ans,sum,tot;
int main()
{
    T=rd();
    register int i;
    while(T--)
    {
        n=rd();mx=0;
        memset(p,0,sizeof(p));
        for(i=1;i<=n;++i) c=rd(),p[c]++,mx=max(mx,c);
        M=(mx+1)<<1;for(m=1,L=-1;m<M;m<<=1) L++;
        for(i=0;i^m;++i) r[i]=(r[i>>1]>>1)|((i&1)<<L);
        for(i=0;i^m;++i) a[i].r=b[i].r=p[i],a[i].i=b[i].i=0;
        fft(a,1);fft(b,1);
        for(i=0;i^m;i++) a[i]=a[i]*b[i];
        fft(a,-1);
        sum=ans=0;tot=(ll)n*(n-1)*(n-2)/6;
        for(i=0;i<=mx;i++)
        {
            sum+=(ll)(a[i].r/m+0.5);
            if(!(i&1)) sum-=p[i>>1];
            if(!p[i]) continue;
            ans+=sum*p[i];
        }
        ans/=2;
        printf("%.7lf
",1.0-ans*1.0/tot);
    }
    return 0;
}