HDU 4609 快速傅里叶变换

题意：给出10^5个边长，问任取三个组成三角形的概率。

这题没遇到还不知道什么是傅里叶变换FFT，当时怎么想时间复杂度都不能降到n*log(n)。实际是这样的，把每个记录长度个数的数组存在num数组中。

直接摘大神的吧 网址http://www.cnblogs.com/kuangbin/archive/2013/07/24/3210565.html

学会了FFT。

这题就很容易了。

其实题目是给了n条线段。问随机取三个，可以组成三角形的概率。

其实就是要求n条线段，选3条组成三角形的选法有多少种。

首先题目给了a数组，

如样例一：

4

1 3 3 4

把这个数组转化成num数组，num[i]表示长度为i的有num[i]条。

样例一就是

num = {0   1   0    2    1}

代表长度0的有0根，长度为1的有1根，长度为2的有0根，长度为3的有两根，长度为4的有1根。

使用FFT解决的问题就是num数组和num数组卷积。

num数组和num数组卷积的解决，其实就是从{1 3 3 4}取一个数，从{1 3 3 4}再取一个数，他们的和每个值各有多少个

例如{0 1 0 2 1}*{0 1 0 2 1} 卷积的结果应该是{0 0  1  0  4  2  4  4  1 }

长度为n的数组和长度为m的数组卷积，结果是长度为n+m-1的数组。

{0 1 0 2 1}*{0 1 0 2 1} 卷积的结果应该是{0 0  1  0  4  2  4  4  1 }。

这个结果的意义如下：

从{1 3 3 4}取一个数，从{1 3 3 4}再取一个数

取两个数和为 2 的取法是一种：1+1

           和为 4 的取法有四种：1+3， 1+3  ，3+1 ，3+1

           和为 5 的取法有两种：1+4 ，4+1；

           和为 6的取法有四种：3+3,3+3,3+3,3+3,3+3

           和为 7 的取法有四种： 3+4,3+4,4+3,4+3

           和为 8 的取法有 一种：4+4

利用FFT可以快速求取循环卷积，具体求解过程不解释了，就是DFT和FFT的基本理论了。

总之FFT就是快速求到了num和num卷积的结果。只要长度满足>=n+m+1.那么就可以用循环卷积得到线性卷积了。

弄完FFT得到一个num数组，这个数组的含义在上面解释过了。

#include <iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
#define N 400005
#define pi acos(-1.0) // PI值
struct complex
{
    double r,i;
    complex(double real=0.0,double image=0.0)
    {
        r=real;
        i=image;
    }
    // 以下为三种虚数运算的定义
    complex operator + (const complex o)
    {
        return complex(r+o.r,i+o.i);
    }
    complex operator - (const complex o)
    {
        return complex(r-o.r,i-o.i);
    }
    complex operator * (const complex o)
    {
        return complex(r*o.r-i*o.i,r*o.i+i*o.r);
    }
} x1[N],x2[N];
void brc(complex *y,int l) // 二进制平摊反转置换 O(logn)
{
    register int i,j,k;
    for(i=1,j=l/2; i<l-1; i++)
    {
        if(i<j)    swap(y[i],y[j]); // 交换互为下标反转的元素
        // i<j保证只交换一次
        k=l/2;
        while(j>=k) // 由最高位检索，遇1变0，遇0变1，跳出
        {
            j-=k;
            k/=2;
        }
        if(j<k)    j+=k;
    }
}
void fft(complex *y,int l,double on) // FFT O(nlogn)
// 其中on==1时为DFT，on==-1为IDFT
{
    register int h,i,j,k;
    complex u,t;
    brc(y,l); // 调用反转置换
    for(h=2; h<=l; h<<=1) // 控制层数
    {
        // 初始化单位复根
        complex wn(cos(on*2*pi/h),sin(on*2*pi/h));
        for(j=0; j<l; j+=h) // 控制起始下标
        {
            complex w(1,0); // 初始化螺旋因子
            for(k=j; k<j+h/2; k++) // 配对
            {
                u=y[k];
                t=w*y[k+h/2];
                y[k]=u+t;
                y[k+h/2]=u-t;
                w=w*wn; // 更新螺旋因子
            } // 据说上面的操作叫蝴蝶操作…
        }
    }
    if(on==-1)    for(i=0; i<l; i++)    y[i].r/=l; // IDFT
}
long long num[N];
long long sum[N];
int a[N>>1];
int main()
{
    int t,n,l,maxn;
    scanf("%d",&t);
    while(t--)
    {
        memset(num,0,sizeof(num));
        maxn=0;
        scanf("%d",&n);
        for(int i=0; i<n; i++)
        {
            scanf("%d",&a[i]);
            num[a[i]]++;
        }
        sort(a,a+n);
        maxn=a[n-1];
        l=1;
        while(l<2*(maxn+1)) l<<=1;
        for(int i=0; i<=maxn; i++)
            x1[i].r=num[i],x1[i].i=0;
        for(int i=maxn+1; i<=l; i++)
            x1[i].r=x1[i].i=0;
        fft(x1,l,1);
        for(int i=0; i<=l; i++) x1[i]=x1[i]*x1[i];
        fft(x1,l,-1);
        for(int i=0; i<=2*maxn; i++) num[i]=(long long)(x1[i].r+0.5);
        for(int i=0; i<n; i++) num[a[i]+a[i]]--;
        for(int i=1; i<=2*maxn; i++) num[i]/=2;
        sum[0]=0;
        for(int i=1; i<=2*maxn; i++)
            sum[i]=sum[i-1]+num[i];
        long long cnt=0;
        for(int i=0; i<n; i++) //枚举每条边并把这条边当成最大边
        {
            cnt+=sum[2*maxn]-sum[a[i]];
            cnt-=n-1;    //减去包括这条边的
            cnt-=(long long)i*(n-i-1);//减去有一条边比这条边大的情况
            cnt-=(long long)(n-i-1)*(n-i-2)/2;//减去有两条都比该条边大的情况
        }
        long long tol=(long long)n*(n-1)*(n-2)/6;
        printf("%.7f\n",double((double)cnt/(double)tol));
    }
    return 0;
}