初探快速傅里叶变换（FFT）

初探快速傅里叶变换，谈谈自己的一些浅见吧，希望以后能把理论基础和相关应用补充一下。

    在信号处理过程中会遇到时域和频域两种不同的表示方法。对于多项式来说 f(x)=a0+a1*x+a2*x^2+....+an*x^n 而言也可以有两种表示方法系数表示法和点集表示法。

   （即相对于对函数而言，给出函数表达式或给出函数经过的顶点两种表达形式）

    （1） 时域——点集表示法

    （2） 频域——系数表示法

    我们通常如果做两个多项式A和B的乘法，当使用系数表示法相乘时与分配率类似，可知复杂度在O(n^2)。

    反观点集表示法(Ax1,Ay1)*（Bx1,By1）可以直接对两个点进行乘法，因此复杂度在O(n)。倘若把多项式的乘法用点集来实现即可降低复杂度。

    快速傅里叶变换（FFT）就是实现在O(nlogn)时间内实现多项式的系数表达和点集表达的相互转化。

    而它的具体实现就是分别选择了(w1,w2.....wn)n个单位复根作为坐标点，根据一些玄学的性质来实现相互转化的（好难待扩充.....）

   

    这样的话当我们遇到子问题包含简单的多项式乘法的问题时，就可以考虑用FFT加速。

    举个栗子好了  (hdu 4609)：题意大概是求给定n个数里任意取三个可以构成三角形的概率。构成三角形需要满足的条件是：两边之和大于第三边。暴力枚举两边显然不行。希望能快速地求出任意两条边的和，我们考虑这样一种构造：设边长为1,2,3,3,5的几条边分别用x,x^2,2*(x^3),x^5来表示，由于多项式的乘法实现的是相同底数时指数的相加，因此可以把构造的一个多项式对自己进行卷积，可以得出一个更长的多项式，其中的次数就分别对应了求出来的和的长度，对应系数则表示有相应的方法数。

    卷积之后得到的结果有什么用呢？对x^k而言，x^0~x^(k-1)的系数和就是所有两边之和小于等于k的方法总数，用前缀和dp处理一下。那么只要排序后枚举第三条最长边i，前面总共i个C(i,2)是总共可能的选边方法，减去小于等于的，就是两边之和大于第三边的合理取法个数。

ps:还是菜到China-Final都没机会去开FFT题....

#include <iostream>
#include <iomanip>
#include <functional>
#include <algorithm>
#include <vector>
#include <string>
#include <queue>
#include <stack>
#include <set>
#include <map>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include <climits>
#include <cctype>
using namespace std;
#define INF 0x3f3f3f3f
#define MP(X,Y) make_pair(X,Y)
#define PB(X) push_back(X)
#define REP(X,N) for(int X=0;X<N;X++)
#define REP2(X,L,R) for(int X=L;X<=R;X++)
#define DEP(X,R,L) for(int X=R;X>L;X--)
#define DEP2(X,R,L) for(int X=R;X>=L;X--)
#define CLR(A,X) memset(A,X,sizeof(A))
#define IT iterator
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> PII;
typedef vector<PII> VII;
typedef vector<int> VI;
#define X first
#define Y second
#define lson(X) ((X)<<1)
#define rson(X) ((X)<<1|1)
const ll maxn=3e5+5;
const double PI=acos(-1.0);
const int M=2e5+5;

struct Complex {
    double r,i;
    Complex(double _r = 0.0,double _i = 0.0) {
        r = _r; i = _i;
    }
    Complex operator +(const Complex &b) {
        return Complex(r+b.r,i+b.i);
    }
    Complex operator -(const Complex &b) {
        return Complex(r-b.r,i-b.i);
    }
    Complex operator *(const Complex &b) {
        return Complex(r*b.r-i*b.i,r*b.i+i*b.r);
    }
};
void change(Complex y[],int len) {
    int i,j,k;
    for(i = 1, j = len/2;i < len-1; i++) {
        if(i < j)swap(y[i],y[j]);
        k = len/2;
        while( j >= k) {
            j -= k;
            k /= 2;
        }
        if(j < k) j += k;
    }
}
 
void FFT(Complex y[],int len,int on) {
    change(y,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 = y[k];
                Complex 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].r /= len;
} 


ll a[maxn],num[maxn],sum[maxn],dp[maxn];
Complex x1[maxn],x2[maxn]; 
int main(){
    int T;
    scanf("%d",&T);
    while(T--){
        CLR(a,0);
        CLR(num,0);
        ll n;
        scanf("%lld",&n);
        ll N=0;
        REP(i,n) {
            scanf("%lld",&a[i]);
            num[a[i]]++;    
            N=max(N,a[i]);
        }
        ll len=1;
        while(len<N*2+1) len<<=1;
        REP2(i,0,N){
            x1[i]=Complex(num[i],0);
            x2[i]=Complex(num[i],0);
        }
        for(int i=N+1;i<len;i++){
            x1[i]=Complex(0,0);
            x2[i]=Complex(0,0);
        }
        FFT(x1,len,1);
        FFT(x2,len,1);
        REP(i,len) x1[i]=x1[i]*x2[i];
        FFT(x1,len,-1);
        REP(i,len){
            sum[i]=(ll)(x1[i].r+0.5);
        //    cout<<sum[i]<<' ';
        }
        //cout<<endl;
        REP(i,n) if(sum[a[i]*2]>0) sum[a[i]*2]--; 
        //REP(i,len) cout<<sum[i]<<' ';cout<<endl;
        dp[0]=0;
        for(ll i=1;i<len;i++) {
            dp[i]=dp[i-1]+sum[i]/2;        
        }
        //REP(i,len) cout<<dp[i]<<' ';cout<<endl;
        //枚举第三条边，C(i,2)-dp[i];
        sort(a,a+n);
        ll ans=0;
        for(ll i=2;i<n;i++){
            ll temp=(ll)(i-1)*(i)/2-dp[a[i]];
            ans+=temp;        
        //    cout<<temp<<' ';
        } 
        //cout<<endl;
        ll tot=n*(n-1)*(n-2)/6;
        //cout<<ans<<' '<<tot<<endl;
        printf("%.7f\n",ans*1.0/tot);
    }
    
    
    
    
    
    return 0;
}
/*
2
4
1 3 3 4
4
2 3 3 4
*/