前言
如果我们能用一种时间上比 \(O(n^2)\) 更优秀的方法来计算大整数（函数）的乘法，那就好了。快速傅里叶变换（FFT） 可以帮我们在 \(O(n\log n)\) 的时间内解决问题。
函数乘积
计算两个大整数之积时，我们发现
23\times45=1035\]
而如果我们把 \((*)\) 式右边的每一位的系数看做一个数每位上的数码，正好得到了 \(1035\)。事实上，对于所有的多项式乘法，以上规律同样成立。
证明： （提示）考虑竖式乘法的过程，和多项式乘法的过程，它们的本质都是一样的。
这样，我们就把问题转换为：计算两个已知函数之积的函数的解析式。
复平面、单位圆
考虑 \(\sqrt{-9}\) 的值。
\]
类似地，\(\forall N\in \Z_-\) 我们都可以用类似的方法得到 $$\sqrt{N}=\sqrt{-N}\times\sqrt{-1}$$
引入虚数单位 \(\text{i}\)，使 \(\text{i}^2=-1.\) 这样我们就重新认识了数的范围，从实数扩充到复数。
一复数 \(a+b\text{i}\) 中的 \(a,b\in\R\)，\(a\) 是它的实数部分，\(b\text{i}\) 是虚数部分。若 \(b=0\)，则它是实数。复数服从实数的大部分运算法则。
若两个复数，它们的实数部分相等，虚数部分之和为 \(0\)，我们称它们互为 共轭复数。
我们知道，数轴上的每个点与每个实数一一对应。类似地，我们可以使用 复平面 上的点表示复数。复平面与平面直角坐标系类似，它的 \(x\) 轴单位长度为 \(1\)，\(y\) 轴单位长度为 \(\text{i}\)。复平面上的点之横纵坐标表示的数之和即为该点表示的数。比如，\((1,2)\) 表示 \(1+2\text{i}\)。
以圆点为圆心，\(1\) 为半径，在复平面上作圆，如图所示。这个圆叫 单位圆。
在圆上任取一点 \(A\)，过此点作 \(x\) 轴的垂线段，垂足为 \(B\)。设 \(∠AOB=\theta\)。易知$$OB=OA\times\cos\theta=\cos\theta\
AB=OB\times\sin\theta=\sin\theta$$\(\therefore A(\cos\theta,\sin\theta).\) 这个 \(\theta\) 叫做 \(A\) 的 辐角。
单位根及其性质
- \(\forall k,n\in\N^*\) 有 \(\omega_n^k\times\omega_n^1=\omega_n^{k+1}\)；
证明：
如图所示，过点 \(A(1,0)\) 作圆内接 \(n\) 边形。设 \(∠AOA'=\theta.\) 则
&=\cos\theta\cos2\theta+(\cos\theta\sin2\theta+\sin\theta\cos2\theta)\times\text{i}-\sin\theta\sin2\theta\\
&=(\sin\theta\sin2\theta+\cos\theta\cos2\theta)+(\sin\theta\cos2\theta+\cos\theta\cos2\theta)\times\text{i}\end{aligned}\]
由三角函数恒等变换式
\cos(\alpha+\beta)=\sin\alpha\cos\beta+\sin\beta\cos\alpha\]
知
&=\sin3\theta+\cos3\theta\times\text{i}\end{aligned}\]
换句话说，\(A'\) 和 \(A''\) 表示的数相乘后得到了对应的 \(A'''\)。
设这个多边形的顶点为 \(\{\omega_i^0\},\ i\in[0,n-1]\)。那么我们一定有$$\omega_n1=\omega_n^{k+1},\ k\in[0,n-1]$$
特别地，我们有 \(\omega_n^n=\omega_n^0\)。它们叫做 \(n\) 次 单位根。
- \(\omega_n^0=\omega_n^n=1\)；
- \(\omega_n^0,\omega_n^1,...,\omega_n^{n-1}\) 互不相同；
- \(\forall k,n\in\N^*\) 有 \(\omega_n^k=\omega_{2n}^{2k}\)；
证明一： 感性理解。一个 \(n\) 边形把单位圆分成了 \(n\) 个部分，取这 \(n\) 个圆弧的中点，顺次连接这 \(2n\) 个点，得到一个新多边形。则 \(\omega_n^k\) 即表示 \(n\) 变形的第 \(k\) 个单位根，也表示 \(2n\) 变形的第 \(2k\) 个点。
证明二：
\(\begin{aligned}\omega_{2n}^{2k}&=\cos2k\times\frac{2\pi}{2n}+\text{i}\times\sin2k\times\frac{2\pi}{2n}\\
&=\cos k\times\frac{2\pi}n+\text{i}\times\sin k\times\frac{2\pi}n\\
&=\omega_n^k.\end{aligned}\)
- \(\forall k,n\in\N^*\) 有 \(\omega_n^{k+\frac2n}=-\omega_n^k\)；
证明一： 感性理解。乘以 \(\omega_n^2\) 的意思其实就是在单位圆逆时针转半圈。单位圆上的一个点，逆时针转了半圈后到达的点，与原来的点关于原点对称。
证明二： $$\begin{aligned}\omega_nk\times\omega_n^{\frac2n}\
&=\omega_n^k\times(\cos\pi+\text{i}\times\sin\pi)\
&=\omega_n^k\times(-1+0)\
&=-\omega_n^k.\end{aligned}$$
- \(\forall k,n\in\N^*\) 有 $$\sum_{i=0}k)^i}=\begin{cases}0,k\neq0\n,k=0\end{cases}$$
证明： \((\text{I})\) 若 \(k\neq0\)，设 $$S=\sum_{i=0}k)^i}\quad...(1)$$ 则 $$\omega_n^k\times S=\sum_{i=1}k)^i}\quad...(2)$$
\((2)-(1)\) 得$$\begin{aligned}(\omega_n^k-1)\times S&=(\omega_nn-1\
S&=\frac{(\omega_nk-1}{\omega_n^k-1}\
&=\frac{1-1}{\omega_n^k-1}\
&=0.\end{aligned}$$
\((\text{II})\) 当 \(k=0\) 时，
\]
## 快速傅里叶变换
显然地，**一个 $n$ 次多项式可以被 $n+1$ 个点唯一确定**。
那么，我们可以在单位圆上取 \(n+1\) 个单位根，代入已知的两个函数，得到 \(n+1\) 对点，再把每对点相乘，得到结果函数上的 \(n+1\) 个点，再求出结果函数。
设已知函数（合并后）为 $$\begin{aligned}f(x)&=\sum_{i=0}i}\
&=(a_0+a_2x+...+a_{n-2}x^{n-2})\
&+(a_1x+a_3x{n-1})\end{aligned}$$
设 $$\begin{aligned}g(x)&=a_0+a_2x+...+a_{n-2}x^{\frac{n}2-1}\
h(x)&=a_1+a_3x+...+a_{n-1}x^{\frac{n}2-1}\end{aligned}$$
则 $$f(x^2)=g(x)+x\times h(x^2)$$
令 \(x=\omega_n^k\)，由 单位根的性质1 得 $$\begin{aligned}f(x{2k})+\omega_n^k\times h(\omega_n^{2k})\
&=g(\omega_{\frac{n}2}k\times h(\omega_{\frac{n}2}^k)\quad...(1)\end{aligned}$$
令 \(x=-\omega_n^k=\omega_n^{k+\frac{n}2}\) 得 $$\begin{aligned}f(x{2k+n})+\omega_n^{k+\frac{n}2}\times h(\omega_n^{2k+n})\
&=g(\omega_n{k+\frac{n}2}\times h(\omega_n^{2k})\
&=g(\omega_{\frac{n}2}k\times h(\omega_{\frac{n}2}^k)\quad...(2)\end{aligned}$$
不难发现，\((1)\) 式与 \((2)\) 式 互轭！换句话说，当你求出 \((1)\) 式的值时，只需将虚数部分取相反数（时间复杂度为 \(O(1)\)）即得到了 \((2)\) 式。
这样，我们就将问题的规模减为一半。同理，剩下的一半也可以使用类似方法，分为更小的两半……没错，这就是 分治，这样就将时间复杂度降为了 \(O(n\log n)\)。
快速傅里叶逆变换
设 \(\{y_i\}\) 是 \(\{a_i\}\) 的傅里叶变换，即在 \(\{\omega_n^i\}\) 处的值；
设 \(\{c_i\}\) 是 \(\{y_i\}\) 在 \(\{\omega_n^{-i}\}\) 的值。
则有 $$\begin{aligned}c_k&=\sum_{i=0}{-k})^i}\
&=\sum_{i=0}{n-1}{a_j(\omega_nj})·(\omega_ni}\
&=\sum_{i=0}{n-1}{(\omega_ni}}\end{aligned}$$
根据 单位根的性质6 ，有且仅有一个 \(j\in[0,n-1]\) 使得 \(j=k\)。
\(\therefore \forall j\) 有且仅有一个$${(\omega_ni}=n$$有且仅有 \(n-1\) 个 $${(\omega_ni}=0$$
&=a_k\times n,\\
a_k&=\frac{c_k}n.\end{aligned}\]
换句话说，经过快速傅里叶逆变换后的数组 \(\{c_i\}\)，除以 \(n\) 后就得到了结果 \(\{a_i\}\)。（这里的 \(n\) 是指 \(c\) 数组的长度）
小结
使用快速傅里叶变换（FFT）计算两个函数之积的步骤如下：
- 分别对两个函数进行快速傅里叶变换；
- 将两组结果相乘；
- 对结果进行快速傅里叶逆变换，并将结果除以 \(n\)。
迭代快速傅里叶变换
经过上文的学习，容易写出以下代码（感谢 linjiayang2016 大佬的代码）
void fast_fast_tle(int limit,complex *a,int type){
if(limit==1)
return ;
complex a1[limit>>1],a2[limit>>1]; //*
for(int i=0;i<=limit;i+=2)
a1[i>>1]=a[i],a2[i>>1]=a[i+1];
fast_fast_tle(limit>>1,a1,type);
fast_fast_tle(limit>>1,a2,type);
complex Wn=complex(std::cos(2.0*Pi/limit),type*std::sin(2.0*Pi/limit)),w=complex(1,0);
for(int i=0;i<(limit>>1);i++,w=w*Wn)
a[i]=a1[i]+w*a2[i],a[i+(limit>>1)]=a1[i]-w*a2[i];
}
发现 \(*\) 行定义了两个不小的数组，而 FFT 函数是被递归调用的，所以会造成爆栈。
| 快速傅里叶变换后序列 | 1 | 5 | 3 | 7 | 2 | 6 | 4 | 8 |
| 原下标的二进制 | 000 | 100 | 010 | 110 | 001 | 101 | 011 | 111 |
| 新下标的二进制 | 000 | 001 | 010 | 011 | 100 | 101 | 110 | 111 |
不难发现，快速傅里叶变换后的序列中，每个数的新下标的二进制，在数值上等于原下标二进制的反转。
设原下标为 \(x\) 的数的新下标为 \(r[x]\)。
以十进制数 \(184\) 的二进制表示为例。设这个二进制数长度为 \(l\)。
\]
我们把这个数分为两部分，第一部分是前 \(l-1\) 位，第二部分是最后一位。不难得到，它的翻转就是 在第二部分后面接上第一部分的翻转形成的数。
即
r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
至此，快速傅里叶变换的相关内容已全部结束。
\(\text{luogu P3803}\)
给你两个正整数 \(n,m\leq10^6\)，和 \(n+1\) 次多项式 \(f(x)\)、\(m+1\) 次多项式 \(g(x)\)。求 \(f(x)\) 和 \(g(x)\) 的卷积。
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
#define reg register
const int MAXN=400010;
const double pi=acos(-1.0);
struct compex{
double x,y;
friend compex operator +(const compex a,const compex b){
compex c;
c.x=a.x+b.x;c.y=a.y+b.y;
return c;
}
friend compex operator -(const compex a,const compex b){
compex c;
c.x=a.x-b.x;c.y=a.y-b.y;
return c;
}
friend compex operator *(const compex a,const compex b){
compex c;
c.x=a.x*b.x-a.y*b.y;
c.y=a.x*b.y+a.y*b.x;
return c;
}
}a[MAXN],b[MAXN];
int r[MAXN];
int n,m,l=0;
int limit=1;
void FFT(compex*t,int type){
for(reg int i=0;i<limit;++i)
if(i<r[i]) swap(t[i],t[r[i]]);
for(reg int mid=1;mid<limit;mid<<=1){
compex wn=(compex){cos(pi/mid),type*sin(pi/mid)};
for(reg int j=0,R=(mid<<1);j<limit;j+=R){
compex w={1.0,0};
for(reg int k=0;k<mid;++k,w=w*wn){
compex x=t[j+k],y=w*t[j+k+mid];
t[j+k]=x+y;
t[j+k+mid]=x-y;
}
}
}
}
int main(){
scanf("%d%d",&n,&m);
for(reg int i=0;i<=n;++i)
scanf("%lf",&a[i].x);
for(reg int i=0;i<=m;++i)
scanf("%lf",&b[i].x);
while(limit<=n+m) limit<<=1,++l;
for(reg int i=1;i<limit;++i)
r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
FFT(a,1);FFT(b,1);
for(reg int i=0;i<limit;++i)
a[i]=a[i]*b[i];
FFT(a,-1);
for(reg int i=0;i<=n+m;++i)
printf("%d ",(int)(a[i].x/limit+0.5));
}