FFT 学习笔记

发布时间 2023-10-31 18:18:09作者: EdGrass

\(FastFuristTransformation\):快速傅立叶变换
——快速求两个多项式的乘积

多项式的点表示法

多项式的性质:用任意\(n+1\)个函数上的不同点均可唯一确定一个多项式。

证明:方程组为一个\(Vandermonder\)行列式,矩阵满秩有唯一解。

截屏2023-10-27 08.07.20
当我们需要多项式 \(A\) 与多项式 \(B\) 的积 \(C\) 的时候,我们只要在 \(A\)\(B\) 中取 \((n+m+1)\) 个点然后直接相乘就可以得到 \(C\) 的点表示法,这个复杂度就是 \(O(n+m+1)\)

复数单位根

将复数域上的单位圆等分成 \(n\) 份,每一份叫做 \(n\)次单位根。
截屏2023-10-30 12.46.37

离散傅立叶正变换

\(DFT\) —— 将系数表示的多项数变为点表示的多项式。

\[ d_k=\sum_{i=0}^{n-1}a_iw_n^{ik} \]

将多项式的奇偶次数分类,变化 \(x\rightarrow x^2\) ,利用复数单位根上的性质作变换后发现,多项式可以由两个 \(\frac{n}{2}\) 长度的多项式得到。然后递归求解的复杂度是 \(O(nlogn)\)

截屏2023-10-30 13.00.46

离散傅立叶逆变换

\(IDFT\) —— 当已知点表示法的一个多项式,得到该多项式系数表示法。

\[ a_k=\frac{1}{n}\sum_{i=0}^{n-1}d_iw_n^{-ki} \]

推导如下:即 \(C_k\)(积多项式的第 \(k\) 项的系数)为 \(\sum{A_iw_n^k{(w_n^{-k})}^i}\)。可得 \(C_k=n*a_k\) ,因为别的项系数都为 \(0\) ,只有 \(k=0\) 这一项系数为 \(n\)
IMG_4106
IMG_4107
然后我们使用 \(FFT\) 优化式子可得:

\[ a_k=\frac{1}{n}\sum_{i=0}^{n}B(w_{n}^{n-k}) \]

这里我们注意到,比较正逆变化的式子只有 \(w_n^k\) 这个 \(k\) 的系数不同,所以几乎是一样的(还有就是多了个系数 \(n\))。

struct comp{
    double r,i;
    comp(double _r=0,double _i=0){r=_r;i=_i;}
    comp operator+ (const comp x){return comp(r+x.r,i+x.i);}
    comp operator- (const comp x){return comp(r-x.r,i-x.i);}
    comp operator* (const comp x){return comp(r*x.r-i*x.i,r*x.i+i*x.r);}
}a[N],b[N];
const double pi=acos(-1.0); 
void FFT(comp a[],int n,int t){
    for(int i=1,j=0;i<n-1;i++){ 
        for(int s=n;j^=s>>=1,~j&s;); 
        if(i<j)swap(a[i],a[j]);
    }
    for(int d=0;(1<<d)<n;d++){
        int m=1<<d,m2=m<<1;
        double o=pi/m*t;
        comp _w(cos(o),sin(o)); 
        for(int i=0;i<n;i+=m2){
            comp w(1,0);
            for(int j=0;j<m;j++){
                comp &A=a[i+j+m],&B=a[i+j],t=w*A;
                A=B-t;B=B+t;w=w*_w; 
            }
        } 
    }
    if(t==-1)for(int i=0;i<n;i++)a[i].r/=n; 
}
int bit=0;
while((1<<bit)<n+m+1)bit++;
bit=(1<<bit);
for(int i=0;i<=n+m;i++){
    cout<<(int)(a[i].r+0.5)<<' ';
}