FFT

原创于 2017-08-21 13:23:27 发布 · 433 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#fft

FFT 专栏收录该内容

6 篇文章

订阅专栏

作者：陈犇

FFT教程

标签（空格分隔）：数学 ACM

FFT是什么？

FFT是快速傅里叶变换(fast Fourier transform)的简称。在ACM领域主要是用来快速求解多项式乘法的算法，
在信号领域也有很大用途

基础知识

卷积
举个例子，给你两个向量 $a (a_0, a_1, a_2), b(b_0, b_1, b_2)$
a和b的卷积就是 $( a_0*b_0, a_1*b_0+a_0*b_1, a_2*b_0+a_1*b_1+a_0*b_2, a_1*b_2+a_2*b_1, a_2*b_2 )$
即可以看作两个多项式 $A(x)=a_0+a_1x^1+a_2x^2和B(x)=b_0+b_1x+b_2x^2$ 相乘。
所以卷积就可看做多项式乘法。那些形如 $c_k = \Sigma_{i=0}^{k}a_ib_{k-i}$ 卷积的就可以类比成多项式乘法，就可通过fft快速求解。
多项式
- 多项式有两种表达形式
- 一种是系数表达形如 $A(x) = \Sigma_{i=0}^{n-1} a_ix^i$
- 另一种是点值表达，形如 $(x_1, y_1), (x_2, y_2), (x_3, y_3)...(x_m, y_m)$ , 且m>=n
- 多项式的点值表示和插值表示可以互推
wn=e2πi/n
- $w_n^{n/2+k}=-w_n^k$
- $w_n^{2k}=w_{n/2}^k$
- $\Sigma_{j=0}^{n-1}w_n^{jk}= (k\%n==0)?n:0$

傅里叶变换的原理

在傅里叶变换里面系数表达和点值表达如何互推？

一个暴力点的方法
系数表达->点值表达 $y_k = \Sigma_{i=0}^{n-1}a_iw_n^{ki}$
点值表达->系数表达 $a_k = 1/n*\Sigma_{i=0}^{n-1}y_iw_n^{-ki}$
证明看黑板
所以不论是系数表达推点值表达还是点值表达推系数表达方法都一样

如何快速的求出上面的式子呢
假设 $A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}$
那 $A(x)=(a_0+a_2x^2+a_4x^4+...+a_{n-2}x^{n-2})+(a_1x+a_3x^3+..+a_{n-1}x^{n-1})$
令 $A_0(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{n/2-1}, A_1(x)=a_1+a_3x+a_5x+...+a_{n-1}x^{n/2-1}$
则 $A(x) = A_0(x^2)+xA_1(x^2)$
那么就可以用 $A_0(w_{n/2}^0), A_0(w_{n/2}^1) ... A_0(w_{n/2}^{n/2-1})$ 和 $A_1(w_{n/2}^0), A_1(w_{n/2}^1) ... A_1(w_{n/2}^{n/2-1}$
得出 $A(w_{n}^0), A(w_{n}^1) ... A(w{n}^{n-1})$

$所以可以采用非递归分值的办法求解FFT$

分治图

void ntt(ll* a, int n, int t) {
    rep(i, 0, n) {
        int rv = rev(i,n);
        if(i<rv) swap(a[i], a[rv]);
    }
    ll g=1;
    if(t==1) g=3;
    else g=inv[3];
    for(int m=2; m<n+1; m<<=1) {
        ll wm= mypow(g, (MOD-1)/m);
        int mid=m>>1;
        for(ll *p=a; p<a+n; p+=m) {
            ll w=1;
            rep(i, 0, mid) {
                ll t=w*(p[mid+i]);
                p[mid+i] = p[i]-t;
                p[i] = p[i]+t;
                w = w*wm%MOD;

                p[i]=%MOD, p[mid+i]%=MOD;
            }
        }
    }
    if(t==-1) {
        rep(i, 0, tn)
            a[i]=a[i]*inv[n]%MOD;
    }
}
void fft(ll *a, ll* b, int n) {
    int tn=MAXN;
    while(tn/4>=n) tn>>=1;
    ntt(a, tn, 1);
    ntt(b, tn, 1);
    rep(i, 0, tn)
        a[i] = a[i]*b[i]%MOD;
    ntt(a, tn, -1);
}