hdu4609 3-idiots（FFT）

题意：

给出n个数枝的长度，
从中随机选出三个，能够组成三角形的概率是多少？

数据范围：n<=1e5，长度a(i)<=1e5

解法：

显然答案=合法方案数/总方案数,
问题变为求合法方案数.

令x^n的系数为长度为n的边的数量,这样就构造出一个多项式,
多项式与自身相乘得到的新多项式就是任意两条边的组合方案数,
要减去自己匹配自己的情况.
同时取法是无序的,所以还要除以2.

因为已知选择两条边长度为n的方案数,O(n)枚举第三条边就能计算出答案.

ll ans=0;//合法方案数
for(int i=1;i<=n;i++){//枚举第三条边a[i],假设a[i]是三条边中最大的,确保只统计一次.
    ans+=sum[len]-sum[[a[i]];//其他两边之和大于a[i]的方案数
    ans-=1ll*(i-1)*(n-i);//减掉x<a[i]<y的情况
    ans-=1ll*(n-i)*(n-i-1)/2;//减掉a[i]<x<y的情况
    ans-=n-1;//减掉x,y中有一个是a[i]的情况,因为a[i]+任何x都>a[i],所以方案数为n-1
}

易错点:
假设a[i]*2的最大值为m,那么是m+1次多项式,所以多项式次数m=a[n]*2+1
ps:
一开始没想到m要+1,然后wa穿了.

code：

#include<bits/stdc++.h>
using namespace std;
const double P=acos(-1.0);
struct CC{//复数
    double x,y;
    CC(double xx=0,double yy=0){x=xx,y=yy;}
    CC operator+(const CC &a)const{return CC(x+a.x,y+a.y);}
    CC operator-(const CC &a)const{return CC(x-a.x,y-a.y);}
    CC operator*(const CC &a)const{return CC(x*a.x-y*a.y,x*a.y+y*a.x);}
};
void change(CC y[],int len){
    for(int i=1,j=len/2;i<len-1;i++){
        if(i<j)swap(y[i],y[j]);
        int k=len/2;
        while(j>=k){
            j-=k;
            k/=2;
        }
        if(j<k)j+=k;
    }
}
void fft(CC y[],int len,int on){//on为1或者-1,-1的时候表示逆变换
    change(y,len);
    for(int h=2;h<=len;h<<=1){
        CC wn(cos(-on*2*P/h),sin(-on*2*P/h));
        for(int j=0;j<len;j+=h){
            CC w(1,0);
            for(int k=j;k<j+h/2;k++){
                CC u=y[k];
                CC t=w*y[k+h/2];
                y[k]=u+t;
                y[k+h/2]=u-t;
                w=w*wn;
            }
        }
    }
    if(on==-1){
        for(int i=0;i<len;i++){
            y[i].x/=len;
        }
    }
}
#define ll long long
const int maxm=6e5+5;
CC x[maxm];
ll cnt[maxm];
ll a[maxm];
ll sum[maxm];
signed main(){
    int T;cin>>T;
    while(T--){
        int n;
        scanf("%d",&n);
        for(int i=1;i<=n;i++){
            scanf("%lld",&a[i]);
        }
        sort(a+1,a+1+n);
        for(int i=1;i<=n;i++){
            cnt[a[i]]++;
        }
        int m=a[n]+1;//len必须是2^k,且len-1>=m,所以这里要加1
        int len=1;
        while(len<m*2)len<<=1;
        for(int i=0;i<=m;i++){
            x[i]=CC(cnt[i],0);
        }
        for(int i=m+1;i<len;i++){
            x[i]=CC(0,0);
        }
        fft(x,len,1);
        for(int i=0;i<len;i++){
            x[i]=x[i]*x[i];
        }
        fft(x,len,-1);
        //
        for(int i=1;i<=len;i++){
            sum[i]=(ll)(x[i].x+0.5);
        }
        len=a[n]*2;
        for(int i=1;i<=n;i++){//自己匹配自己的情况
            sum[a[i]*2]--;
        }
        for(int i=1;i<=len;i++){//无序,除以2
            sum[i]/=2;
        }
        sum[0]=0;
        for(int i=1;i<=len;i++){
            sum[i]+=sum[i-1];
        }
        //
        ll tot=1ll*n*(n-1)*(n-2)/6;//总方案数C(n,3)
        ll ans=0;//合法方案数
        for(int i=1;i<=n;i++){//枚举第三条边a[i],假设a[i]是三条边中最大的,确保只统计一次.
            ans+=sum[len]-sum[a[i]];//其他两边之和大于a[i]的方案数
            ans-=1ll*(i-1)*(n-i);//减掉x<a[i]<y的情况,一个取左边一个取右边
            ans-=1ll*(n-i)*(n-i-1)/2;//减掉a[i]<x<y的情况,两个都取右边
            ans-=n-1;//减掉x,y中已经有一个是a[i]的情况,因为a[i]+任何x都>a[i],所以方案数为n-1
        }
        printf("%.7f\n",1.0*ans/tot);
        //
        for(int i=1;i<=n;i++){//clear
            cnt[a[i]]=0;
        }
    }
    return 0;
}

---