FFT学习记录（雾）

前置知识：

1.点值表示

设$A(x)$是一个$(n-1)$次多项式，那么把$(n)$个不同的$(x)$代入，会得到$(n)$个$(y)$。这$(n)$对$(x,y)$唯一确定了该多项式，即只有一个多项式能同时满足“代入这些$(x)$，得到的分别是这些$(y)$”。

由多项式可以求出其点值表示，而由点值表示也可以求出多项式。

2.复数知识

复数相乘的规则：模长相乘，幅角相加。

3.单位根的性质

性质一：${\omega }_{2n}^{2k}={\omega }_n^k$

证明：它们对应的点/向量是相同的。

性质二：${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$

证明：它们对应的点是关于原点对称的（对应的向量是等大反向的）。

1.什么是FFT

快速傅里叶变换（FFT）是一种能在$O(n\log n)$的时间内将一个多项式转换成它的点值表示的算法

2.为什么要使用FFT

FFT可以用来加速多项式乘法

假设有两个$n-1$次的多项式$A(x)$和$B(x)$，我们的目标是把它们乘起来。

普通的多项式乘法是$O(n^2)$的。

因为我们需要枚举$A(x)$中的每一项，分别与$B(x)$中的每一项相乘得到多项式$C(x)$

但有趣的是，两个用点值表示的多项式相乘，复杂度是$O(n)$的。

即：

$C(x_i)=A(x_i)\times B(x_i)$，所以$O(n)$枚举$x_i$就行了

如果把两个多项式转换成点值表示，再相乘，在把新的点值表示转换成多项式不就可以$O(n)$解决多项式乘法了？

但是，把多项式转换成点值表示的朴素算法是$O(n^2)$。

所以大整数乘法的复杂度瓶颈就在于多项式转换成点值表示这一步，只要完成这步，就能$O(n)$求答案了。

这时傅里叶就跑出来了…

3.离散傅里叶变换（快速傅里叶变换的朴素版）

傅里叶发明了一种办法：规定点值表示中的$n$个$x$为$n$个模长为$1$的复数。

C++的STL提供了复数模板

定义：complex<double> x;
运算：直接使用加减乘除。

傅里叶要用到的$n$个复数，不是随机找的，而是——把单位圆（圆心为原点、$1$为半径的圆）$n$等分，取这$n$个点（或点表示的向量）所表示的虚数，即分别以这$n$个点的横坐标为实部、纵坐标为虚部，所构成的虚数。

从点$(1,0)$开始（该点显然是要去的点之一)，逆时针将这$n$个点从$0$开始编号，第$k$个点对应的虚数记作${\omega }_n^k$（根据复数相乘时模长相乘幅角相加可以看出，${\omega }_n^k$是${\omega }_n^1$的$k$次方，所以${\omega }_n^1$被称为$n$次单位根）。

根据每个复数的幅角，可以计算出所对应的点/向量。${\omega }_n^k$对应的点/向量是$(\cos \frac{k}{n}2\pi ,\sin \frac{k}{n}2\pi )$，也就是说这个复数是$\cos \frac{k}{n}2\pi +i\sin \frac{k}{n}2\pi$

傅里叶说：把$n$个复数${\omega }_n^0,{\omega }_n^1,{\omega }_n^2,\dots ,{\omega }_n^{n-1}$代入多项式，能得到一种特殊的点值表示，这种点值表示就叫离散傅里叶变换吧！

为什么要使用单位根作为$x$代入

设$(y_0,y_1,y_2,\dots ,y_{n-1})$为多项式$A(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_{n-1}{x}^{n-1}$的离散傅里叶变换。

再设一个多项式$B(x)=y_0+y_{1}x+y_2{x}^2+\cdots +y_{n-1}{x}^{n-1}$

现在我们把上面的$n$个单位根的倒数，即
${\omega }_n^0,{\omega }_n^{-1},{\omega }_n^{-2},\dots ,{\omega }_n^{-(n-1)}$作为$x$代入到$B(x)$，得到了一个新的离散傅里叶变换$(z_0,z_1,z_2,\dots ,z_{n-1})$。

$$\begin{array}{rl}{z}_k& =\sum _{i=0}^{n-1}{y}_i{\left({\omega }_n^{-k}\right)}^i\\ & =\sum _{i=0}^{n-1}\left(\sum _{j=0}^{n-1}{a}_j{\left({\omega }_n^i\right)}^j\right){\left({\omega }_n^{-k}\right)}^i\\ & =\sum _{j=0}^{n-1}{a}_j\left(\sum _{i=0}^{n-1}{\left({\omega }_n^{j-k}\right)}^i\right)\end{array}$$

其中${\sum }_{i=0}^{n-1}\left({\omega }_n^{j-k}\right)$是可求的：当$j-k=0$时，它等于$n$。其他时候通过等比数列和可知它等于
$$\frac{{\left({\omega }_n^{j-k}\right)}^n-1}{{\omega }_n^{j-k}-1}=\frac{{\left({\omega }_n^n\right)}^{j-k}-1}{{\omega }_n^{j-k}-1}=\frac{1^{j-k}-1}{{\omega }_n^{j-k}-1}=0$$
那么$z_k$就等于$n{a}_k$，即：
$$a_i=\frac{z_i}{n}$$

一个结论

把多项式$A(x)$的离散傅里叶变换结果作为另一个多项式$B(x)$的系数，取单位根的倒数即${\omega }_n^0,{\omega }_n^{-1},{\omega }_n^{-2},...,{\omega }_n^{-(n-1)}$作为$x$代入$B(x)$，得到的每个数再除以$n$，得到的是$A(x)$的各项系数。

这实现了傅里叶变换的逆变换——把点值表示转换成多项式系数表示，这就是离散傅里叶变换神奇的性质。

4.快速傅里叶变换

虽然傅里叶发明了神奇的变换，能把多项式转换成点值表示又转换回来，但是……它仍然是暴力代入的做法，复杂度仍是$O(n^2)$

于是，快速傅里叶变换应运而生。它是一种分治的傅里叶变换。

快速傅里叶变换的数学证明

同样的，设$A(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_{n-1}{x}^{n-1}$，现在为了求离散傅里叶变换，要把一个$x={\omega }_n^k$代入。

考虑将$A(x)$的每一项按照下标奇偶分成两部分：

$A(x)=\left(a_o+a_2{x}^2+\cdots +a_{n-2}{x}^{n-2}\right)+\left(a_{1}x+a_3{x}^3+\cdots +a_{n-1}{x}^{n-1}\right)$

设两个多项式：
$$A_1(x)=a_0+a_{2}x+...+a_{n-2}{x}^{\frac{n}{2}-1}$$

$$A_2(x)=a_1+a_{3}x+...+a_{n-1}{x}^{\frac{n}{2}-1}$$

则：
$$A(x)=A_1(x^2)+x{A}_2(x^2)$$
假设$k<\frac{n}{2}$，现在将$x={\omega }_n^k$代入：
$$\begin{array}{rl}A\left({\omega }_n^k\right)& =A_1\left({\omega }_n^{2k}\right)+{\omega }_n^k{A}_2\left({\omega }_n^{2k}\right)\\ & =A_1\left({\omega }_{\frac{n}{2}}^k\right)+{\omega }_n^k{A}_2\left({\omega }_{\frac{n}{2}}^k\right)\end{array}$$
那么对于$A\left({\omega }_n^{k+\frac{n}{2}}\right)$
$$\begin{array}{rl}A\left({\omega }_n^{k+\frac{n}{2}}\right)& =A_1\left({\omega }_n^{2k+n}\right)+{\omega }_n^{k+\frac{n}{2}}{A}_2\left({\omega }_n^{2k+n}\right)\\ & =A_1\left({\omega }_{\frac{n}{2}}^k\times {\omega }_n^n\right)+{\omega }_n^{k+\frac{n}{2}}{A}_2\left({\omega }_{\frac{n}{2}}^k\times {\omega }_n^n\right)\\ & =A_1\left({\omega }_{\frac{n}{2}}^k\right)-{\omega }_n^k{A}_2\left({\omega }_{\frac{n}{2}}^k\right)\end{array}$$
所以，如果我们能知道两个多项式$A_1(x)$和$A_2(x)$分别在$\left({\omega }_{\frac{n}{2}}^0,{\omega }_{\frac{n}{2}}^1,{\omega }_{\frac{n}{2}}^2,\dots ,{\omega }_{\frac{n}{2}}^{\frac{n}{2}-1}\right)$的点值表示，就可以$O(n)$求出$A(x)$在${\omega }_n^0,{\omega }_n^1,{\omega }_n^2,\cdots {\omega }_n^{n-1}$处的点值表示了。而$A_1(x)$和$A_2(x)$都是规模缩小了一半的子问题。分治边界是$n=1$，此时直接return。

快速傅里叶变换的实现

typedef complex<double> cp;
const double PI=acos(-1);
cp omega(int n,int k)
{
    return cp(cos(2*PI*k/n),sin(2*PI*k/n));
}
void fft(cp  *a,int n,bool inv)
{
    if(n==1) 
    {
        return ;
    }
    static cp buf[N];
    int m=n/2;
    for(int i=0;i<m;i++)//将每一项按照奇偶分为两组
    { 
        buf[i]=a[2*i];
        buf[i+m]=a[2*i+1];
    }
    for(int i=0;i<n;i++)
    {
        a[i]=buf[i];
    }
    fft(a,m,inv); //递归处理两个子问题
    fft(a+m,m,inv);
    for(int i=0;i<m;i++)//枚举x，计算A(x)
    { 
        cp x=omega(n,i); 
        if(inv) 
        {
            x=conj(x);//conj是一个自带的求共轭复数的函数，精度较高。当复数模为1时，共轭复数等于倒数
        }
        buf[i]=a[i]+x*a[i+m];//根据之前推出的结论计算
        buf[i+m]=a[i]-x*a[i+m];
    }
    for(int i=0;i<n;i++)
    {
        a[i]=buf[i];
    }
}

inv表示这次用的单位根是否要取倒数。