【文文殿下】快速傅里叶变换(FFT)学习笔记

多项式

定义

形如(A(x)=sum_{i=0}^{n-1} a_i x^i)的式子称为多项式。
我们把(n)称为该多项式的次数界。
显然，一个(n-1)次多项式的次数界为(n)。

运算法则

设(A(x))和(B(x))为多项式，且次数界分别为(n),(m)。则有：
(A(x)=sum_{i=0}^{n-1}a_i x^i)
(B(x)=sum_{i=0}^{m-1}b_i x^i)

他们遵循下面的常用运算法则：

(A(x)+B(x)=sum_{i=0}^{max(n-1,m-1)} (a_i+b_i)x^i)
(A(x)-B(x)=sum_{i=0}^{max(n-1,m-1)} (a_i-b_i)x^i)
(A(x)B(x)=sum_{k=0}^{n+m-2}sum_{i=0}^{k}a_i b_{k-i} x^k)

两个多项式的乘积称为他们的卷积。卷积也是一个多项式。
易证次数界分别为(n),(m)的多项式的卷积是一个次数界为(n+m-1)的多项式。
用定义计算卷积的复杂度为(O(n^2))
使用快速傅里叶变换，可以使复杂度降为(O(nlogn))

复数

代数形式

形如(a+bi)形式的数称为复数。其中(i)是(sqrt{-1})
我们可以使用平面上的点((a,b))来描述一个复数。
复数的运算法则和向量相似。

三角形式

除了上文提到的代数形式，也有一种三角形式来表示复数。
如图所示。

我们知道复数(a+bi)对应着平面上的点((a,b))，也对应着复平面上的一个向量。
这个向量的长度叫做复数(a+bi)的模，记作(|a+bi|)，通常用小写字母(r)表示。
这个向量与(x)正半轴有一个夹角，称为辐角，我们用希腊字母( heta)表示。

显然有(a=rcos heta),(b=rsin heta)
把上式带入到代数表达式，可以得到(a+bi=rcos heta+irsin heta=r(cos heta+isin heta))

定理:两复数相乘，模长相乘，辐角相加。
证明：
设(z_1=r_1(cos x+isin x)),(z_2=r_2(cos y+isin y)).
则(z_1 z_2=r_1 r_2 [cos xcos y-sin xsin y+i(cos xsin y+cos ysin x)]=r_1 r_2[cos(x+y)+isin(x+y)])

指数形式

欧拉公式：(e^{i heta}=cos heta+isin heta)
证明：
将(e^x)，(sin x)，(cos x)按照泰勒展开得到他们的泰勒级数：
(e^x=1+frac{x}{1!}+frac{x^2}{2!}+frac{x^3}{3!}+…)
(cos x=1-frac{x^2}{2!}+frac{x^4}{4!}-frac{x^6}{6!}+…)
(sin x=x-frac{x^3}{3!}+frac{x^5}{5!}-frac{x^7}{7!})
将(x=i heta)带入得：
(e^{i heta}=1+ heta i-frac{ heta^2}{2!}-frac{ heta^3}{3!}i+frac{ heta^4}{4!}+frac{ heta^5}{5!}i-frac{ heta^6}{6!}-frac{ heta^7}{7!}i+…)
(e^{i heta}=(1-frac{ heta^2}{2!}+frac{ heta^4}{4!}-frac{ heta^6}{6!}+…)+i( heta-frac{ heta^3}{3!}+frac{ heta^5}{5!}-frac{ heta^7}{7!}+…))
即(e^{i heta}=cos heta+isin heta)
将欧拉公式带入到三角形式的表达式中就能得到复数的指数形式表达式:(z=r(cos heta+isin heta)=re^{i heta})

单位复数根

在复数域内，有且仅有(n)个数，使得这个数的(n)次方等于1，我们把这样的数称为单位复数根，记作(omega_{n}^{i}　(iepsilon [0,n-1]))
把( heta=2pi)带入欧拉公式，得(e^{2pi i}=1)
所以(omega_{n}^{k}=(e^{frac{2pi i}{n}})^k)
消去引理：(omega_{dn}^{dk}=omega_{n}^{k})
证明：(omega_{dn}^{dk}=(e^{frac{2pi i}{dn}})^{dk}=(e^{frac{2pi i}{n}})^k=omega_{n}^{k})

折半引理：((omega_{n}^{k})^2=omega_{n/2}^{k})
证明：将(d=frac {1}{2})带入消去引理

求和引理：(sum_{i=0}^{n-1} (omega_{n}^{k})^i=0) ((k)不是(n)的倍数)
证明：利用等比数列求和公式，原式等于(frac{(omega_{n}^{k})^n-1}{omega_{n}^{k}-1}=frac{0}{omega_{n}^{k}-1})
因为分母不能为0，所以(omega_{n}^{k})不等于(1),所以(k)不是(n)的倍数.

循环性质:((omega_{n}^{k})^p=omega_{n}^{pk}=omega_{n}^{(pkmod{n})})
证明：(omega_{n}^{kp}=omega_{n}^{(kpmod n+left lfloor frac{kp}{n} ight floor n)}=omega_{n}^{(kp mod n)})

对称性质：(w_{n}^{k+frac{n}{2}}=-w_{n}^{k})
证明：(w_{n}^{k+frac{n}{2}}=e^{frac{2pi ki}{n}+pi i}=w_{n}^{k} e^{pi i}=-w_{n}^{k})

离散傅里叶变换

多项式的两种表示方法

开头提到的多项式，都是用系数表达法来表示的。
也就是将次数界为(n)的多项式的系数(a_0,a_1,…,a_{n-1})看成(n)维向量(vec{a}=(a_0,a_1,…,a_{n-1})),(vec{a})称作多项式(A(x))的系数表示。
事实上，除了系数表示(A(x)=sum_{i=0}^{n-1}a_i x^i)还有一种表示方法：点值表示法。
一个次数界为(n)的多项式，就是(n)的点值对组成的集合。((x_0,y_0),(x_1,y_1),…,(x_{n-1},y_{n-1})) .对于(0<=k<=n-1)，所有的(x_k)各不相同，记作(y_k=A(x_k))
而离散傅里叶变换就是要快速求出多项式在这(n)的点的值

Cooley-Tukey算法

Cooley-Tukey算法是一种基于分治的算法
们取(n)个(n)次单位根 (w_{n}^{0},w_{n}^{1},…,w_{n}^{n-1}) 进行求值:(A(w_{n}^{k})=sum_{i=0}^{n-1}a_i w_{n}^{ki})
向量(vec{y}=(A(w_{n}^{0}),…,A(w_{n}^{n-1})))称作(vec{a}=(a_0,…,a_n-1))的离散傅里叶变换
接下来我们按照奇偶分开来算：
(A(w_{n}^{k})=sum_{i=0}^{n-1}a_i w_{n}^{ki})
(=sum_{i=0}^{frac{n}{2}-1}a_{2i}w_{n}^{2ki}+w_{n}^{k}sum_{i=0}^{frac{n}{2}-1}a_{2i+1}w_{n}^{2ki})
(=sum_{i=0}^{frac{n}{2}-1}a_{2i}w_{frac{n}{2}}^{ki}+w_{n}^{k}sum_{i=0}^{frac{n}{2}-1}a_{2i+1}w_{frac{n}{2}}^{ki}) (折半引理)
而(A(w_{n}^{k+frac{n}{2}})=sum_{i=0}^{frac{n}{2}-1}a_{2i}w_{frac{n}{2}}^{ki}-w_{n}^{k}sum_{i=0}^{frac{n}{2}-1}a_{2i+1}w_{frac{n}{2}}^{ki})(对称性质)
这样问题就变成了求(frac{n}{2})个单位根的点值，问题的规模缩小了一半，所以是(O(nlogn))的。

快速傅里叶变换的逆变换

求值与插值

对于一个多项式，如果知道它的系数表达式，我们很容易求出它的点值表达式，这一过程成为求值
如果知道它的点值表达式，让你求它的系数表达式，这一过程称为插值
定理：当插值多项式的次数界等于已知点值对的数目时，插值才是明确的
证明：(y_k=A(x_k))等价于

其中，左边的矩阵称为范德蒙矩阵，其行列式(Det=prod_{0leq j<kleq n-1}(x_k-x_j))
如果所有(x_k)均不相等，那么其行列式不为(0)，该矩阵可逆。
因此给定点值表达，能够唯一确定系数表达
所以如果我们知道了一个多项式的点值表示，就可以通过矩阵求逆的方法来算出系数向量
现在我们得到了点值表达式，那么如何求系数表达式呢？
矩阵求逆的复杂度是(O(n^3))的,范德蒙矩阵求逆可以优化至(O(n^2))，但是复杂度依旧很高.
我们考虑构造一个现成的能算的东西进行优化

构造逆矩阵

我们现在已经得到了这样的方程：

记系数矩阵为(V)，记下面的矩阵(d_{ij}=w_{n}^{-ij})

它们相乘的结果(E=DV)满足(e_{ij}=sum_{k=0}^{n-1}w_{n}^{k(j-i)})
当(i=j)时，(e_{ij}=n)
当(i eq j)时，(e_{ij}=sum_{k=0}^{n-1}(w_{n}^{(j-i)})^k=0)(求和引理)
所以(frac{1}{n}D)就是我们要求的逆矩阵，如下图：

这就相当于我们把单位复数根(w^{k}_{n})换成(w^{-k}_{n})再做一遍离散傅里叶变换结果除以(n)就行了

模板代码

#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
const int maxn = 1200000;
const double pi = acos(-1.0);
struct Complex {
    double r,i;
    Complex(double r,double i):r(r),i(i){}
    Complex(){};
};
Complex operator + (Complex a,Complex b) {
    return Complex(a.r+b.r,a.i+b.i);
}
Complex operator - (Complex a,Complex b) {
    return Complex(a.r-b.r,a.i-b.i);
}
Complex operator * (Complex a,Complex b) {
    return Complex(a.r*b.r-a.i*b.i,a.r*b.i+b.r*a.i);
}
int n,m,L[maxn],R[maxn];
Complex A[maxn<<1],B[maxn<<1];
inline void read(int &x) {
    char ch=getchar();x=0;
    while(ch<'0'||ch>'9')ch=getchar();
    while(ch>='0'&&ch<='9') {x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
}
inline void read(double &x) {
    char ch=getchar();x=0;
    while(ch<'0'||ch>'9')ch=getchar();
    while(ch>='0'&&ch<='9') {x=10*x+ch-'0';ch=getchar();}
}
void fft(Complex *a,int n,int inv) {
    for(register int i=1,j=n/2;i<n-1;++i) {
        if(i<j) std::swap(a[i],a[j]);
        int k = n>>1;
        while(j>=k) {
            j-=k;k>>=1;
        }
        j+=k;
    }
    for(register int j=2;j<=n;j<<=1) {
        Complex wn(cos(2*pi/j*inv),sin(2*pi/j*inv));
        for(register int i=0;i<n;i+=j) {
            Complex w(1,0);
            for(register int k=i;k<i+(j>>1);k++) {
                Complex u=a[k],t=a[k+(j>>1)]*w;
                a[k]=u+t;a[k+(j>>1)]=u-t;
                w=w*wn;
            }
        }
    }
    if(inv==-1) for(register int i=0;i<n;i++) a[i].r=a[i].r/n+0.5;
}
int ans[maxn<<1];
int main() {
    read(n);read(m);
    for(register int i=0;i<=n;i++) read(A[i].r);
    for(register int i=0;i<=m;i++) read(B[i].r);
    int limit = 1;
    while(limit<=n+m) limit<<=1;
    
    fft(A,limit,1);fft(B,limit,1);

    for(register int i=0;i<=limit;i++) {
        A[i]=A[i]*B[i];
    }
    fft(A,limit,-1);
    for(register int i=0;i<=m+n;i++) printf("%llu ",(unsigned long long int)(A[i].r));
    return 0;
}