快速傅里叶变换（FFT）学习笔记

多项式

有 $n$ 次多项式可以形式化的写为：

$f(x)=\sum\limits_{i=0}^{n}a_ix^i$

其中序列 $a$ 可以叫做这个多项式的系数序列。

在下文，我们会用 $f_i$ 表示多项式 $f$ 的系数序列 $a$ 的第 $i$ 项，即为 $a_i$。

多项式的表示方法

• 系数表示法

显然，只要我们得知了系数序列，就可以唯一确定一个多项式。

• 点值表示法

对于一个 $n$ 次多项式，代入 $n+1$ 个互不相同的值 $x_i$，可以得到 $n+1$ 个满足 $y_i=f(x_i)$ 的 $y_i$。这时，序列 $x$ 和序列 $y$ 可以唯一确定 $f$。

因为有了 $n+1$ 个点 $(x_i,y_i)$，就可以通过插值来确定 $f$。

多项式的运算

设 $f,g$ 是两个 $n$ 次多项式，有：

$(f+g)(x)=\sum\limits_{i=0}^n(f_i+g_i)x^i$

$(f*g)(x)=\sum\limits_{k=0}^{2n}x^k\sum\limits_{i\le n,j\le n,i+j=k}f_ig_j$

多项式除法可以使用多项式求逆做。

容易发现，在系数表示法下：

• 多项式加法是 $O(n)$ 的，因为直接系数相加就行

• 多项式乘法是 $O(n^2)$ 的，因为需要枚举 $i,j$

而在点值表示法下：（注意此时的运算需要满足序列 $x$ 相同）

• 多项式加法是 $O(n)$ 的，因为直接相加就行（乘法分配律）

• 多项式乘法是 $O(n)$ 的，因为直接让 $y$ 相乘就行

所以在点值表示法下做乘法运算是很优的！

但是直接通过代入把系数表示法转换为点值表示法是 $O(n^2)$ 的，而 FFT 就是把这一过程优化到了 $O(n\log n)$。

$n$ 次本原单位根

FFT 能优化时间复杂度的一个重要原因是它代入的不是随随便便的 $n+1$ 个不同的值，而是 $n$ 次本原单位根的 $0$ 次方、$1$ 次方一直到 $n$ 次方。

$n$ 次单位根其实就是任意一个就是满足 $x^n=1$ 的 $x$，这个东西在实数范围显然只有不多于两个，但是在虚数范围就有 $n$ 个了。

而 $n$ 次本原单位根则是一个特殊的 $n$ 次单位根 $\omega_n$，满足 $\omega_n^0\not=\omega_n^1\not=\omega_n^2\not=\dots\not=\omega_n^n$。

想要构造一个 $n$ 次本原单位根，最好的方法就是把复平面上的单位圆平分成 $n$ 份。例如六次本原单位根 $\omega_6$：

显然，$n$ 次本原单位根 $\omega_n=\operatorname{cos}(\dfrac{2\pi}{n})+\operatorname{sin}(\dfrac{2\pi}{n})i$。

我们考虑 $\omega_n^2$ 是什么。因为虚数的乘法法则是模相乘，辐角相加，而 $\omega_n$ 的模是 $1$，所以 $\omega_n^2$ 相当于 $\omega^n$ 转了一下：

所以 $\omega_n^k=\omega_n^{k\mod n}$，相当于转一圈后会转回来。本原单位根会转是一个很重要的性质！（虽然暂时没有用 ）

DFT 和 IDFT

现在我们可以把单位根们代入多项式求值了，这就是 DFT。设 $f$ 是一个 $n-1$ 次多项式，那么有：

$\hat f_k=\sum\limits_{i=0}^{n-1}f_i\cdot\omega_n^{ki}$

下面让我们证明 DFT 的逆变换 IDFT：

$f_k=\dfrac{1}{n}\sum\limits_{i=0}^{n-1}\hat f_i\cdot\omega_n^{-ki}$

将 DFT 的柿子代入：

$f_k=\dfrac{1}{n}\sum\limits_{i=0}^{n-1}\omega_n^{-ki}\sum\limits_{j=0}^{n-1}f_j\cdot\omega_n^{ij}$

$=\dfrac{1}{n}\sum\limits_{j=0}^{n-1}f_j\sum\limits_{i=0}^{n-1}\omega_n^{-ki}\cdot\omega_n^{ij}$

$=\dfrac{1}{n}\sum\limits_{j=0}^{n-1}f_j\sum\limits_{i=0}^{n-1}\omega_n^{i(j-k)}$

考虑 $\sum\limits_{i=0}^{n-1}\omega_n^{i(j-k)}$ 这部分：

• 若 $j=k$，那么 $\sum\limits_{i=0}^{n-1}\omega_n^{i(j-k)}=\sum\limits_{i=0}^{n-1}\omega_n^{0}=n$

• 若 $j\not=k$，那么因为 $0\le j,k<n$，所以 $|j-k|<n,\omega_n^{j-k}\not=0$，由等比数列求和公式得到：

$\sum\limits_{i=0}^{n-1}\omega_n^{i(j-k)}=\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i$

$=\dfrac{1-(\omega_n^{j-k})^n}{1-\omega_n^{j-k}}$

$=\dfrac{1-(\omega_n^n)^{j-k}}{1-\omega_n^{j-k}}$

$=\dfrac{1-1}{1-\omega_n^{j-k}}$

$=0$

所以

$f_k=\dfrac{1}{n}\sum\limits_{j=0}^{n-1}f_j\cdot n\cdot[j=k]$

$=\dfrac{1}{n}\cdot n\cdot f_k$

$=f_k$

得证。

可以发现 DFT 和 IDFT 的柿子惊人地相似，只不过是多了个 $\dfrac{1}{n}$ 和本原单位根指数上面的负号而已。所以我们只需要解决 DFT，IDFT 就能迎刃而解。

FFT

设 $f$ 是一个次数为 $n-1$ 且 $n$ 是偶数的多项式（如果原来的 $n$ 是奇数那么可以补一项系数为 $0$ 的项），设 $m=\dfrac{n}{2}$，那么有：

$f(x)=\sum\limits_{i=0}^{n-1}f_ix^i$

根据奇偶性分类，有：

$f(x)=\sum\limits_{i=0}^{m-1}f_{2i}x^{2i}+\sum\limits_{i=0}^{m-1}f_{2i+1}x^{2i+1}$

$=\sum\limits_{i=0}^{m-1}f_{2i}x^{2i}+x\sum\limits_{i=0}^{m-1}f_{2i+1}x^{2i}$

定义 $f0(x)=\sum\limits_{i=0}^{m-1}f_{2i}x^{i}$，$f1(x)=\sum\limits_{i=0}^{m-1}f_{2i+1}x^{i}$，那么有：

$f(x)=f0(x^2)+x \cdot f1(x^2)$

在代入单位根之前，我们先来看两个单位根的性质：

$\omega_{2n}^{2k}=\omega_{n}^k$

$\omega_{2n}^{n+k}=-\omega_{2n}^k$

第一条性质是因为以 $\omega_{2n}$ 的辐角转 $2k$ 次相当于以 $\omega_{n}$ 的辐角转 $k$ 次。

而第二条性质是因为以 $\omega_{2n}$ 的辐角转 $n$ 次就到达了 $-1$ 的位置。

记住这两条性质，代入单位根，对于满足 $0\le k<m$ 的 $k$，有：

$f(\omega_n^k)=f0(\omega_n^{2k})+\omega_n^k\cdot f1(\omega_n^{2k})$

$=f0(\omega_m^k)+\omega_n^k\cdot f1(\omega_m^k)$

$f(\omega_n^{m+k})=f0(\omega_n^{2(m+k)})+\omega_n^{m+k}\cdot f1(\omega_n^{2(m+k)})$

$=f0(\omega_m^k)-\omega_n^{k}\cdot f1(\omega_m^k)$

以上两个柿子十分相似，所以被称作蝴蝶操作。

由蝴蝶操作不难得出，只要我们算出了 $f0(\omega_m^0,\omega_m^1,\omega_m^2,\dots,\omega_m^{m-1})$ 和 $f1(\omega_m^0,\omega_m^1,\omega_m^2,\dots,\omega_m^{m-1})$，就可以快速得到 $f(\omega_n^0,\omega_n^1,\omega_n^2,\dots,\omega_n^{n-1})$。

而 $f0(\omega_m^0,\omega_m^1,\omega_m^2,\dots,\omega_m^{m-1})$ 和 $f1(\omega_m^0,\omega_m^1,\omega_m^2,\dots,\omega_m^{m-1})$ 又可以递归下去算，所以我们可以分治下去，用 $O(n\log n)$ 的时间复杂度完成 DFT！

不过需要注意的是，$n$ 必须是 $2$ 的幂，要不然无法分治。

快速进行 IDFT 也相当好办，只要把 FFT 中用到的所有 $\omega_n$ 换成 $\omega_n^{-1}$ ，最后再把所有点值乘上 $\dfrac{1}{n}$ 即可。

但是这还不够快，所以会被卡常。观察到影响效率的是按奇偶性分治，所以可以考虑令 $\operatorname{rev}(x)$ 表示 $x$ 的二进制反过来组成的数，令 $f^\prime_i=f_{\operatorname{rev}(i)}$，那么在 $f^\prime$ 上的分治操作其实就是相邻两个序列的分治，这样可以大大加快 FFT 的速度，并且可以把 FFT 写成非递归的。

代码如下：（模板题）

#include <iostream>
#include <cstdio>
#include <cmath>

using namespace std;

const long long MS=5000005;
const double PI=acos(-1);

struct plex
{
    double x,y;
    plex(double a=0,double b=0) {x=a,y=b;}
};
plex operator+(plex a,plex b) {return plex(a.x+b.x,a.y+b.y);}
plex operator-(plex a,plex b) {return plex(a.x-b.x,a.y-b.y);}
plex operator*(plex a,plex b) {return plex(a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y);}

inline int getlen(int n)
{
    int res=1;
    while(res<n) res<<=1;
    return res;
}
int p_rev[MS],p_rev_lstn;
inline void FFT(int n,plex a[],int tpe)
{
    if(p_rev_lstn!=n)
    {
        p_rev_lstn=n;
        for(int i=0;i<n;i++) p_rev[i]=(p_rev[i>>1]>>1)|((i&1)?n>>1:0);
    }
    for(int i=0;i<n;i++) if(p_rev[i]<i) swap(a[p_rev[i]],a[i]);
    for(int mid=1;mid<n;mid<<=1)
    {
        int len=mid<<1;
        plex Wn=plex(cos(2*PI/len),tpe*sin(2*PI/len));
        for(int l=0;l<n-len+1;l+=len)
        {
            plex Wk=plex(1,0);
            for(int k=0;k<mid;k++,Wk=Wk*Wn)
            {
                plex x=a[l+k],y=Wk*a[l+mid+k];
                a[l+k]=x+y,a[l+mid+k]=x-y;
            }
        }
    }
}
inline void DFT(int n,plex a[]) {FFT(n,a,1);}
inline void IDFT(int n,plex a[])
{
    FFT(n,a,-1);
    for(int i=0;i<n;i++) a[i].x/=n,a[i].y/=n;
}

int n,m;
plex a[MS],b[MS];

int main()
{
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++)
	{
		scanf("%lf",&a[i].x);
	}
	for(int i=0;i<=m;i++)
	{
		scanf("%lf",&b[i].x);
	}
	int len=getlen(n+m+1);
	DFT(len,a);
	DFT(len,b);
	for(int i=0;i<len;i++)
	{
		a[i]=a[i]*b[i];
	}
	IDFT(len,a);
	for(int i=0;i<n+m+1;i++)
	{
		printf("%d ",(int)(a[i].x+0.5));
	}
	printf("\n");
	return 0;
}