2013 多校第一场 hdu 4609 3-idiots（FFT）

hdu 4609

题目：http://acm.hdu.edu.cn/showproblem.php?pid=4609

题目大意：给你 n 条边，问你任意三条边能组成三角形的概率。

思路：设 Ai 记录长度为 i 的边的个数，让它对自己本身做FFT，求出任意两条边能得到的和的数量。然后枚举第三条边，两边之和大于第三条边，第三条边指的是三条边中最长的边，所以还要 - （比它大的两条边的组合+一大一小+含有它自身的组合）。

以上是看了解题报告后写的，自己先开始只能想到时间复杂度为O（n^2）的DP，明显超时。借着这道题学了下 FFT ，看了一个下午加一个晚上才终于懂了。建议大家学的话，最好看算法导论，其他书写的都不是很详细，网上大部分也写的很乱。代码很多也是参考别人的代码写的，然后又自己按照思路重新敲了一遍 = =

纪念自己的第一道传说中的 F！F！T！

代码如下：

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace  std;

const double PI = acos(-1.0);

typedef __int64 lld;

const int MAXN = 111111 ;

struct Com
{
    double re,im;
    Com(double a =0,double b=0)
    {
        re = a;
        im = b;
    }
    Com operator + (const Com &tmp)
    {
        return Com(re+tmp.re,im+tmp.im);
    }
    Com operator - (const Com &tmp)
    {
        return Com(re-tmp.re,im-tmp.im);
    }
    Com operator * (const Com &tmp)
    {
        return Com(re*tmp.re-im*tmp.im,re*tmp.im+im*tmp.re);
    }
};

void change(Com y[] , int len)//雷德算法
{
    for(int i = 1,j = len>>1;i<len-1;i++)
    {
        if(i<j) swap(y[i],y[j]);
        int k = len>>1;
        while(j>=k)
        {
            j -= k;
            k >>= 1;
        }
        j+=k;
    }
}

void fft(Com y[],int len,int on)//on == 1 表示正变换，-1 表示反变换
{
    change(y,len);
    for(int h = 2;h <= len;h<<=1)
    {
        int m = h;
        Com wn = Com(cos(on*2*PI/m),sin(on*2*PI/m));
        for(int i = 0;i<len;i+=m)
        {
            Com w = Com(1,0);
            for(int j = 0;j<m/2;j++)
            {
                Com u = y[i+j];
                Com t = w*y[i+j+m/2];
                y[i+j] = u+t;
                y[i+j+m/2] = u-t;
                w = w*wn;
            }
        }
    }
    if(on == -1)
    {
        for(int i=0;i<len;i++)
            y[i].re /= len;
    }
}

int a[MAXN];
Com x1[MAXN<<2];//最坏情况可能会接近于4倍
lld num[MAXN<<2],sum[MAXN<<2];//10^5*10^5 会爆int

int main()
{
    int T;
    scanf("%d",&T);
    while(T--)
    {
        int n;
        scanf("%d",&n);
        memset(num,0,sizeof(num));
        for(int i=0;i<n;i++)
        {
            scanf("%d",&a[i]);
            num[a[i]]++;
        }
        sort(a,a+n);
        int len_a = a[n-1]+1;
        int len = 1;
        while(len<2*len_a) len<<=1;//2n次多项式,长度为2n
        for(int i=0;i<len_a;i++)
            x1[i] = Com(num[i],0);
        for(int i = len_a;i<len;i++)
            x1[i] = Com(0,0);//补0

        fft(x1,len,1);
        for(int i = 0;i<len;i++)
            x1[i] = x1[i]*x1[i];
        fft(x1,len,-1);

        for(int i = 0;i<len;i++)
            num[i] = (lld)(x1[i].re+0.5);
        for(int i = 0;i<n;i++)
            num[a[i]+a[i]]--;//减去自己与自己的组合
        for(int i = 0;i<len;i++)
            num[i] /= 2;//考虑a+b,b+a的组合，个数/2
        sum[0]=0;
        for(int i = 1;i<len;i++)
            sum[i] = sum[i-1]+num[i];//求前项和
        lld cnt = 0;
        for(int i = 0;i < n; i++)
        {
            cnt += sum[len-1]-sum[a[i]];

            cnt -= (n-1);//减去自身和其他线段的组合
            cnt -= (lld)i*(n-i-1);//减去一大一小的组合
            cnt -= (lld)(n-i-1)*(n-i-2)/2;//减去两个都比它大的组合
        }
        //printf("%I64d\n",cnt);
        lld tot = (lld)n*(n-1)*(n-2)/6;
        printf("%.7f\n",(double)cnt/tot);
    }
    return 0;
}