FFT学习笔记

FFT可用于解决一些卷积问题。
一般问题形式如下：
$C=A*B$
$C[i]=\sum_{j=0}^iA[i]*B[i-j]$

若把A，B看成两个次数为n多项式
$A(x)=\sum_{i=0}^na[i]*x^i,B(x)=\sum_{i=0}^nb[i]*x^i$
原问题等于两个多项式相乘，C的次数等于2n-1

点值

一个次数界为n的多项式A的点值表达为n个点值对所组成的集合。形如
{(x0,y0),(x1,y1)…(xn−1,yn−1)}
其中yi=A(xi)，x各不相同。

插值

插值运算为点值运算的逆运算，如果给定一个n个点值对的点值表达，我们可以确定唯一一个次数界为n的多项式。

对于多项式乘法，
我们可以求出A的点值表达{(x0,a0),(x1,a1)…(xn−1,an−1)}和B的点值表达{(x0,b0),(x1,b1)…(xn−1,bn−1)}
则C的点值表达为{(x0,a0*b0),(x1,a1*b1)…(xn−1,an−1*bn-1)}
再通过插值运算，则可求出C的每一项的系数。

直接暴力，复杂度是O(n^2)的，FFT的关键就在于x的选取。

复数

形如a+bi($i=\sqrt{-1}$)
如果建立平面直角坐标系，a作为x坐标，b作为y坐标，我们可以用一个点对(a,b)来表示一个复数，也相当于一个平面向量。
用三角函数简单推导，可以发现，复数乘法得到的新向量的模长等于原向量的模长之积，与x轴的夹角等于原向量的夹角之和。

n次单位复数根

n次单位复数根是满足$w^n=1$的复数w
且$w^1$~$w^n$把平面均匀的分成n份
记$w$为主n次单位复数根，则w=(cos($2\pi/n$),sin($2\pi/n$))
$w^i=w^{i-1}*w=(cos(2i\pi/n$),sin($2i\pi/n))$
为了区分其他次单位复数根，我们添一个下标来表示($w^i_n$)

消去引理

$w^{dk}_{dn}=w^k_n$
$w^{dk}_{dn}=(cos(2dk\pi/(dn)),sin(2dk\pi/(dn)))=(cos(2k\pi/n),sin(2k\pi/n))=w^k_n$

折半引理

2n次单位复数根的平方的集合，等于n次单位复数根的集合。
易知对于任意一个<=n的k，存在唯一一个<=n的k’满足$w^k_n=-w^{k'}_n$
$(w^k_{2n})^2=(w^{k'}_{2n})^2=w^k_n$

求和引理

k!=0，n>0时$\sum_{i=0}^{n-1}w^{ki}_n=0$
相当于等比数列求和，$\sum_{i=0}^{n-1}w^{ki}_n=(w^{kn}_n-w^0_n)/(w^k-1)=0$

FFT

为了方便，我们只考虑次数n=2^k的多项式乘法
如何求出单个数x的函数值A(x)？我们定义两个多项式
$A0(x)=a0+a2x+a4x^2⋅⋅⋅an−2x^{n/2}$
$A1(x)=a1+a3x+a5x^2⋅⋅⋅an−1x^{n/2}$
$A(x)=A0(x^2)+x*A1(x^2)$
原问题：求A(x)在n次单位复数根上的函数值。
转化成：求A0(x)和A1(x)在n/2次单位复数根上的值。
于是可以递归地对n/2的多项式A0(x)与A1(x)在n/2个n/2次单位复数根进行求值。

插值运算
如果把点值运算写成矩阵方程的形式，可以得到表达式$Y=V_n*A$。
那么插值运算相当于求A，$A=V^{-1}_n*Y$

结论，$V^{−1}_n[j,k]=w^{−jk}_n/n$.
证明：$[V_n∗V^{−1}_n][j,j′]=\sum_{i=0}^{n-1}w^{(j′−j)i}_n/n$，如果j’=j，那后面为1，否则根据求和引理，为0。

实现优化

如何把空间复杂度降为O(n)呢？
不好讲，详见代码吧…

struct z{
    double x,y;
    z(double _x=0,double _y=0){x=_x,y=_y;}
}c[maxn],d[maxn],g[maxn],a[maxn],b[maxn],an[maxn];
z operator *(z a,z b){return z(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
z operator +(z a,z b){return z(a.x+b.x,a.y+b.y);}
z operator -(z a,z b){return z(a.x-b.x,a.y-b.y);}
void dft(ar *a,int sig){
    int i,m=1,k,half,j;
    fo(i,0,l-1) g[f[i]]=a[i];
    while (m<l){
        half=m,m+=m;
        fo(k,0,half-1){
            ar w=ar(cos(pi*sig*k/half),sin(pi*sig*k/half));
            for(i=k;i<l;i+=m){
                j=i+half;
                ar u=g[i],v=w*g[j];
                g[i]=u+v,g[j]=u-v;
            }
        }
    }fo(i,0,l-1)a[i]=g[i];
    if (sig<0) fo(i,0,l-1) a[i].x/=l;
}
void fft(){
    fo(i,0,l-1) a[i]=b[i]=z(0,0);
    fo(i,0,n) a[i].x=a1[i],b[i].x=b1[i];
    dft(a,1),dft(b,1);
    fo(i,0,l-1) an[i]=a[i]*b[i];
    dft(an,-1);
    fo(i,0,l-1) ans[i]=(ll)(an[i].x+eps);
}
void chu(){
    l=1,t=0;
    while (l<=n*2) l+=l,t++;
    fo(i,0,l-1){
        int x=i,y=t-1;f[i]=0;
        while (x) f[i]+=(x&1)<<y,y--,x>>=1;
    }
}