HDU 4609 3-idiots

题意：

给出n条边，问选择3个边能构成三角形的概率

题解：

边权不大

考虑FFT的组合意义，得到的f[i]的值，是所有的j+k=i的a[j]*b[k]

如果把下标当做边权的话，FFT一次，就得到了选择两个边，能组成i的方案数。

三角形充要条件，最小的两个边之和大于最大边。

但是最大边可能之前在较小的用过了。

所以，考虑不合法的情况，最大边c>=a+b

这样，c一定就是一个>=a+b的后缀和了。

C(n,3)减去不合法的即可。

对于任意一个合法的，一定没有被减去，任意一个不合法的一定减去了。

正确。

注意：

1.多组数据清空，，用到的多项式是m>2*mx，buc的上界是mx

2.buc*buc可能爆int

3.FFT注意精度误差，绝对值可能在(0~0.5)之间，不能直接取整，采用四舍五入

a=floor(b+0.5)或者，a=round(b)

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cstring>
#include<iostream>
#include<cmath>
#define reg register int
#define il inline
#define numb (ch^'0')
using namespace std;
typedef long long ll;
il void rd(int &x){
    char ch;x=0;bool fl=false;
    while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true);
    for(x=numb;isdigit(ch=getchar());x=x*10+numb);
    (fl==true)&&(x=-x);
}
namespace Miracle{
const int N=1e5+5;
const double Pi=acos(-1);
ll n;
ll buc[4*N],sum[4*N];
struct node{
    double x,y;
    node(){}
    node(double xx,double yy){
        x=xx;y=yy;
    }
    node friend operator *(node a,node b){
        return node(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
    }
    node friend operator +(node a,node b){
        return node(a.x+b.x,a.y+b.y);
    }
    node friend operator -(node a,node b){
        return node(a.x-b.x,a.y-b.y);
    }
    void clear(){
        x=0.0,y=0.0;
    }
}a[4*N],b[4*N],c[4*N];

int rev[4*N];
void fft(node *f,int n,int c){
    for(reg i=0;i<n;++i){
        if(i<rev[i]) swap(f[i],f[rev[i]]);
    }
    for(reg p=2;p<=n;p<<=1){
        int len=p/2;
        node gen=node(cos(Pi/len),c*sin(Pi/len));
        for(reg l=0;l<n;l+=p){
            node buf=node(1,0);
            for(reg k=l;k<l+len;++k){
                node tmp=buf*f[k+len];
                f[k+len]=f[k]-tmp;
                f[k]=f[k]+tmp;
                buf=buf*gen;
            }
        }
    }
}
ll mx;
int m;
void clear(){
    
    for(reg i=0;i<m;++i) a[i].clear(),b[i].clear(),c[i].clear();
    for(reg i=0;i<=mx;++i) buc[i]=0;mx=0;
}
int main(){
    int T;
    rd(T);
    while(T--){
        scanf("%lld",&n);
        clear();
        int x;
//        for(reg i=0;i<16;++i){
//            cout<<buc[i]<<" ";
//        }cout<<endl;
        for(reg i=1;i<=n;++i){
            rd(x);mx=max(mx,(ll)x);
            buc[x]++;
        }
        m=1;
        for(;m<=mx*2;m<<=1);
        for(reg i=0;i<m;++i){
            rev[i]=(rev[i>>1]>>1)|((i&1)?(m>>1):0);
            a[i].x=buc[i];
            b[i].x=buc[i];
        }
        //cout<<" mx "<<mx<<" "<<m<<endl;
        //cout<<" buc "<<endl;
//        for(reg i=0;i<m;++i){
//            cout<<buc[i]<<" ";
//        }cout<<endl;
        fft(a,m,1);fft(b,m,1);
        for(reg i=0;i<m;++i) b[i]=a[i]*b[i];
        fft(b,m,-1);
        //cout<<" dughrg "<<endl;
        for(reg i=0;i<m;++i) b[i].x=round(b[i].x/m);
        sum[mx]=buc[mx];
        for(reg i=mx-1;i>=0;--i) sum[i]=sum[i+1]+buc[i];
        double tot=n*(n-1)*(n-2)/6;
        double ans=tot;
        for(reg i=1;i<=mx;++i){
            //cout<<" ii "<<i<<" "<<b[i].x<<endl;
            ll tmp=0;
            if(i%2==0){
                tmp=buc[i/2];
            }
            b[i].x-=tmp*tmp;
            b[i].x/=2;
            b[i].x+=tmp*(tmp-1)/2;
            //cout<<" b[i] "<<b[i].x<<endl;
            ans-=b[i].x*sum[i];    
        }
        printf("%.7lf
",ans/tot);
    }
    return 0;
}

}
signed main(){
    Miracle::main();
    return 0;
}

/*
   Author: *Miracle*
   Date: 2018/12/21 21:03:28
*/