FFT详解

多项式乘法:

$$C[i]=A[j]*B[i-j]$$
可以在O(n*m)时间内求出

多项式点值表示

n次多项式可以用n+1个二维平面上的点(x,y)来唯一表示

证明

设f(n)为一个n次多项式,g(n)为另一个1次多项式，且其点值表示相同
设其点值表示为x[n+1]
那么有f(x[i])-g(x[i])==0
令l(i)为f(i)-g(i)
可知l(i)零点有n+1个,与代数基本定理（一个 n-1 次多项式在复数域上有且仅有 n-1 个根）相矛盾矛盾
所以点值表示可以唯一表示一个多项式
注:下文的n并不是实际的n,而是>=实际n的第一个形如2^k的数，前面不够的项数可以补0,不影响多项式乘法的计算

DFT离散傅里叶变换

->将多项式的系数表示变为点值表示
直接计算虽然有快速的秦九韶算法，但要求n+1个数的对应函数值，时间复杂度为O(n^2)

FFT快速傅里叶变换

复数

形如x+bi的数为复数,x为其实数部分,bi为其虚数部分,i代表$\sqrt{-1}$

复平面

x轴表示复数的实部,y轴表示复数的虚部的二维平面,称作复平面

复数可以用复平面上的向量表示上图的向量就表示复数3+i

单位根

单位根，表示为 $e^{i\theta}$ ，可知 $e^{i\theta}=cos\theta+i\cdot sin\theta$

用$\omega_n^1$表示n的一次单位根,也就是满足将单位圆n等分的最小正角度对应的向量
也就是$e^{2\pi\over n}$
$\omega_n^n==1$

性质一

$\omega_n^i=\omega_n^j*\omega_n^{i-j}$
根据定义$\omega_n^i=e^{2\pi i\over n}$,还是很显然的

性质二

$\omega_{2n}^{2i}=\omega_n^i$
根据定义$2\pi 2i\over 2n$,上下同时消去2

性质三

$\omega_{n}^{i}=-\omega_n^{i+{n\over2}}$
由$\omega_{n}^{i+{n\over2}}=cos{2\pi (i+{n\over2})\over{ n}}+i\cdot sin{2\pi (i+{n\over2})\over{ n}}$
$=cos({2\pi i\over{ n}}+\pi)+i\cdot sin({2\pi i\over{ n}}+\pi)$
$=-cos{2\pi i\over{ n}}-i*sin{2\pi i\over{ n}}$
$=-\omega_{n}^{i}$

FFT的实现

考虑将单位根次方带入多项式求点值
对于多项式A(x)
可以将其奇数项和偶数项分开
有

$$A(x)==a_0+a_2*x^2+...a_n*x^n+x*(a_1+a_3*x^2...+a_{n-1}*x^{n-2})$$
令奇数项部分为 $A1(x)=a_1+a_3*x^..+a_{n-1}*x^{{n-2}\over2}$,
偶数项为 $A2(x)=a_0+a_2*x^2+...a_n*x^{n\over2}$
那么有 $$A(x)==A2(x^2)+x*A1(x^2)$$
带入单位根里就是
$$A(\omega_{n}^{i})=A2(\omega_{n\over2}^{i})+\omega_{n}^{i}*A1(\omega_{n\over2}^{i})$$
考虑 $i<n/2$的情况 $$A(\omega_{n}^{i})=A2(\omega_{n\over2}^{i})+\omega_{n}^{i}*A1(\omega_{n\over2}^{i})$$
$$A(\omega_{n}^{i+{n\over2}})=A2(\omega_{n\over2}^{i+{n\over2}})+\omega_{n}^{i+{n\over2}}*A1(\omega_{n\over2}^{i+{n\over2}})$$
$$A(\omega_{n}^{i+{n\over2}})=A2(\omega_{n\over2}^{i+{n\over2}})+\omega_{n}^{i+{n\over2}}*A1(\omega_{n\over2}^{i+{n\over2}})$$
$$A(\omega_{n}^{i+{n\over2}})=A2(\omega_{n}^{2i+{n}})-\omega_{n}^{i}*A1(\omega_{n}^{2i+{n}})$$
$$A(\omega_{n}^{i+{n\over2}})=A2(\omega_{n}^{2i})-\omega_{n}^{i}*A1(\omega_{n}^{2i})$$
$$A(\omega_{n}^{i+{n\over2}})=A2(\omega_{n\over2}^{i})-\omega_{n}^{i}*A1(\omega_{n\over2}^{i})$$
所以只要知道前一半的数,就可以对应推出后一半的数
可以采用倍增算法

IDFT离散傅里叶反变换

将点值表示重新计算成系数表示
假设di为DFT后$\omega_{n}^{i}$对应的点值
对di建立多项式,将$\omega_{n}^{-k}$代入,算出点值表示
令$C_k$为$\omega_{n}^{-k}$对应的点值表示
则$C_k=\sum_{i=0}^{n-1}d_i*(\omega_{n}^{-k})^i$
$C_k=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j*(\omega_{n}^{i})^j*(\omega_{n}^{-k})^i$
$C_k=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j*(\omega_{n}^{i})^{j-k}$
$C_k=\sum_{j=0}^{n-1}a_j*\sum_{i=0}^{n-1}(\omega_{n}^{i})^{j-k}$
令j-k=l
令$s(j,k)=\sum_{i=0}^{n-1}(\omega_{n}^{i})^{j-k}$
若$j==k$
$s(j,k)==n$
若$j!=k$
$s(j,k)=\omega_{n}^{0}+\omega_{n}^{(j-k)}+\omega_{n}^{2*(j-k)}+...+\omega_{n}^{(n-1)*(j-k)}$
是公比为$\omega_{n}^{(j-k)}$的等比数列
即$s(j,k)=\omega_{n}^{0}*{{{\omega_{n}^{(j-k)}}^{n}-1}\over{\omega_{n}^{(j-k)}-1}}$
$s(j,k)={{{\omega_{n}^{n}}^{(j-k)}-1}\over{\omega_{n}^{(j-k)}-1}}$
$s(j,k)={{1^{(j-k)}-1}\over{\omega_{n}^{(j-k)}-1}}=0$
所以$C_k=\sum_{j=0}^{n-1}a_j*[j==k]*n$
即$C_k=a_k*n$
即$C_k/n=a_k$
所以可以算出$C_k$再倒推$a_k$

实现

考虑递归实现
只需要在这一层从左儿子区间推出这个区间内容
需要O(n)的辅助空间
但常数比较大
考虑从下面层向上递推

其实是按照二进制位反过来排序，再按堆结构递推即可

总结

处理多项式乘法$C(x)=A(x)*B(x)$时
先对A(x)求一下点值表达（DFT），再对B(x)求一下点值表达(DFT)，由于点值是可以直接相乘的，可以O(n)做乘法
再对最后求得的C(x)的点值表达做一遍IDFT求出系数表达
时间复杂度O(3*nlogn+n*2)
也就是O(nlogn)

代码：

inline void FFT(C *a,int opt){
    for(int i=0;i<N;i++)
        if(i<r[i])
            swap(a[i],a[r[i]]);//二进制位排序
    for(int i=1;i<N;i<<=1){//子区间长度
        C wn;
        wn.a=cos(opt*pii/i),wn.b=sin(opt*pii/i);
        for(int j=0;j<N;j+=(i<<1)){//区间起点
            C w,x,y;
            w.a=1,w.b=0;
            for(int k=0;k<i;k++,w=w*wn){
                x=a[j+k],y=w*a[i+j+k];
                a[j+k]=x+y;
                a[j+k+i]=x-y;
            }
        }
    }
    if(opt==-1)//IDFT
        for(int i=0;i<N;i++){
            a[i].a/=N;
        }
}

例题:
COGS1473. 很强的乘法问题
设A[i]为一个乘数从后往前数第i位上的数(从0开始)
B[i]为另一个乘数,C[i]为积
由于乘法满足C[i+j]=A[i]*B[j]
对每一位上的数看成系数,做多项式乘法
最后在O(n)进位
代码：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <cmath>
using namespace std;
const int maxn=(1<<20);
const double pii=acos(-1);
struct C{
    double a,b;
    C(){a=0,b=0;}
    C operator *(const C &x)const{
        C y;
        y.a=a*x.a-b*x.b;
        y.b=a*x.b+b*x.a;
        return y;
    }
    C operator +(const C &x)const{
        C y;
        y.a=x.a+a;
        y.b=x.b+b;
        return y;
    }
    C operator -(const C &x)const{
        C y;
        y.a=a-x.a;
        y.b=b-x.b;
        return y;
    }
};
C A[maxn],B[maxn];
int r[maxn];
int ans[maxn];
string x,y;
int n,m;
int N;
inline void FFT(C *a,int opt){
    for(int i=0;i<N;i++)
        if(i<r[i])
            swap(a[i],a[r[i]]);
    for(int i=1;i<N;i<<=1){
        C wn;
        wn.a=cos(opt*pii/i),wn.b=sin(opt*pii/i);
        for(int j=0;j<N;j+=(i<<1)){
            C w,x,y;
            w.a=1,w.b=0;
            for(int k=0;k<i;k++,w=w*wn){
                x=a[j+k],y=w*a[i+j+k];
                a[j+k]=x+y;
                a[j+k+i]=x-y;
            }
        }
    }
    if(opt==-1)
        for(int i=0;i<N;i++){
            //printf("%d\n",int(a[i].a));
            a[i].a/=N;
        }
}
int main(){
    freopen("bettermul.in","r",stdin);
    freopen("bettermul.out","w",stdout);
    cin>>x>>y;
    n=x.length(),m=y.length();
    for(int i=0;i<n;i++){
        A[i].a=x[n-i-1]-'0';
        //printf("%d\n",int(A[i].a));
    }
    for(int i=0;i<m;i++){
        B[i].a=y[m-i-1]-'0';
        //printf("%lf\n",B[i].a);
    }
    n=n+m-1;
    int L=0;
    while((1<<L)<n)
        L++;
    N=(1<<L);
    //printf("%d\n",N);
    for(int i=0;i<N;i++){
        r[i]=(r[i>>1]>>1)|((i&1)?(N>>1):0);
        //printf("%d\n",r[i]);
    }
    FFT(A,1);
    FFT(B,1);
    for(int i=0;i<N;i++)
        A[i]=A[i]*B[i];
    FFT(A,-1);
    for(int i=0;i<n;i++)
        ans[i]=int(A[i].a+0.5);
    for(int i=0;i<n;i++){
        ans[i+1]+=ans[i]/10;
        ans[i]%=10;
        if(ans[n]&&i==n-1)
            n++;
    }
    for(int i=n-1;i>=0;i--)
        printf("%d",ans[i]);
return 0;
}

感谢

邓祎明的zhihu