HDU-4609 3-idiots（FFT）

转载请注明出处，谢谢http://blog.csdn.net/bigtiao097?viewmode=contents
参考文章http://blog.csdn.net/jackyguo1992/article/details/12622139

题意：

给n条边，问随便取三个可以组成三角形的概率

思路：

先求出任意取两条不同的边其和为i的情况数，即为数组num[i]，这个需要用到FFT（这里不解释FFT的原理），因为1e5的数据范围，$n^2$ 复杂度会超时,FFT可以在$nlog(n)$的时间内求出， 用FFT求出来的数组需要减去取同一个的情况，并且要除以2（取a、b和取b、a是同一种情况）
求出总的方案数$tot = \frac {n\times(n-1)\times(n-2)}{6}$，注意这里一定是选取了三条不一样的边，这样在后面减去多余的情况是只需要减去任意两边之和不能大于第三边的情况就好了，接下来就是枚举第三条边，减去$a+b<=c$（三角形任意两边之和大于第三边） 的情况

---

具体代码如下：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const double PI = acos(-1);
const int maxn = 1e5+5;//maxn = max(len1,len2)
const int maxm = 4*maxn;//k是>=maxn的最小的2的幂,maxm=2*k
int xx[maxn], yy[maxn], len1, len2,len;//len1，len2分别为两个多项式的最高次数+1
int num[maxm];
int a[maxn];
int n;
//复数结构体
struct Complex
{
    double x, y;//实部为x，虚部为y
    Complex(double x=0, double y=0):x(x),y(y){}
    Complex operator+(const Complex &rhs) const { return Complex(x+rhs.x, y+rhs.y);}
    Complex operator-(const Complex &rhs) const { return Complex(x-rhs.x, y-rhs.y);}
    Complex operator*(const Complex &rhs) const { return Complex(x*rhs.x-y*rhs.y,x*rhs.y+y*rhs.x);}
};
Complex x1[maxm],x2[maxm];
/*
进行FFT和IFFT前的反转变换。
将位置i和(i二进制反转后位置)互换。
len必须取2的幂
*/
void change(Complex *x, int len)
{
    Complex t;
    for(int i = 1, j = len/2; i < len-1; i++)
    {
        //交换下标互反的元素，i<j保证交换一次
        //i做正常的+1，j做反转的+1，始终保持i和j是反转的
        if(i < j)
        {
            t = x[i];
            x[i] = x[j];
            x[j] = t;
        }
        int k = len / 2;
        while(j >= k)
        {
            j -= k;
            k >>= 1;
        }
        if(j < k) j += k;
    }
}
/*
做FFT
len必须为2的幂
on==1时是DFT，on==-1时是IDFT
*/
void fft(Complex *x, int len, int on)
{
    change(x, len);
    for(int h = 2; h <= len; h <<= 1)
    {
        Complex wn(cos(-on*2*PI/h), sin(-on*2*PI/h));
        for(int j = 0; j < len; j += h)
        {
            Complex w(1, 0);
            for(int k = j; k < j+h/2; k++)
            {
                Complex u = x[k];
                Complex t = w*x[k+h/2];
                x[k] = u+t;
                x[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if(on == -1)
        for(int i = 0; i < len; i++) x[i].x /= len;
}
void workFFT()
{
    len = 1;
    memset(x1,0,sizeof x1);
    memset(x2,0,sizeof x2);
    memset(num,0,sizeof num);
    while(len < len1*2 || len < len2*2) len <<=1;
    for(int i = 0; i < len; i++)
            x1[i] =  x2[i] = Complex(0,0);
    for(int i = 0; i < len1; i++) x1[i] = Complex(xx[i],0);
    for(int i = 0; i < len2; i++) x2[i] = Complex(yy[i],0);
    fft(x1,len,1);
    fft(x2,len,1);
    for(int i = 0; i < len; i++) x1[i] = x1[i]*x2[i];
    fft(x1,len,-1);
    for(int i = 0; i < len; i++) num[i] = (int)(x1[i].x+0.5);
}
int main()
{
    int T;
    cin>>T;
    while(T--)
    {
        scanf("%d",&n);
        len1 = 0;
        for(int i=1;i<=n;i++)
            scanf("%d",&a[i]);
        memset(xx,0,sizeof xx);
        memset(yy,0,sizeof yy);
        for(int i=1;i<=n;i++)
        {
            xx[a[i]]++;
            yy[a[i]]++;
        }
        sort(a+1,a+1+n);
        len1 = a[n]+1;
        len2 = len1;
        workFFT();
        for(int i=1;i<=n;i++)
            num[a[i]+a[i]]--;
        for(int i=1;i<=len;i++)
            num[i]/=2;

        ll tot = 1LL*n*(n-1)*(n-2)/6;
        for(int i=1,j=0;i<=len;i++)
            if(num[i])
        {
            while(j<=n&&a[j]<i)j++;
            tot-=1LL*num[i]*(n-j+1);
        }
        printf("%.7f\n",tot*1.0*6/n/(n-1)/(n-2));
    }
}