快速傅里叶变换FFT
快速傅里叶变换FFT(Fast Fourier Transformation)
本文主要讲述如何使用FFT来实现快速多项式乘法。
多项式的表示
系数表示
对于一个多项式
向量\(a=(a_0, a_1, ..., a_{n-1})\)为多项式的系数表示。
利用系数表示时，给出\(x_0\)，在\(O(n)\)的时间内可以求出\(A(x_0)\)的值，而且\(A(x)\)的加减法也可以在\(O(n)\)内求出，但乘法则需要\(O(n^2)\)的时间。
点值表示
一个多项式\(A(x)\)的点值表示为一个由\(n\)个点值对所组成的集合
使得对\(k=0, 1, 2, ..., n-1\), 所有\(x_k\)各不相同，且\(y_k=A(x_k)\).一个多项式可以有很多不同的点值表示。
次数界
如果一个多项式\(A(x)\)的最高次的非零系数为\(a_k\)，则称\(A(x)\)的次数为\(k\)，记\(degree(A)=k\)。任何严格大于一个多项式次数的整数都是该多项式的次数界。
定理：对于任意\(n\)个点值对组成的集合，其中\(x_k\)都不同，那么存在唯一的次数界为\(n\)的多项式\(A(x)\)，满足\(y_k=A(x_k), k=0, 1, 2, ..., n-1\)。
证明：由线性代数中的范德蒙行列式可证。
点值表示的优点
相对于系数表示，点值表示的加减乘法都可以在\(O(n)\)内算出(其中乘法因为当\(x\)一定时，\(y\)值相乘就是之后的\(x\)对应的\(y\)值。
因为一个多项式的点值表示可以有多种，所以可以选取一些较容易算的值，然后快速把一个多项式从系数表示转化为点值表示，运用点值表示相乘，然后再快速从点值表示转化为系数表示。
DFT
因此选择单位复数根作为\(x\)，其中\(n\)次单位复数根是满足\(\omega ^n=1\)的复数\(\omega\)，\(n\)次单位复数根有\(n\)个：\(\omega_k=e^{2\pi ik/n}, k=0, 1, 2, ..., n-1, \omega_n=e^{2\pi i/n}\)称为主\(n\)次单位根。
引理1（消去引理）：对于任何整数\(n, k \geq 0, d>0, \omega_{dn}^{dk}=\omega_{n}^{k}\)，特别地，对于任意偶数\(n>0\)，有\(\omega_{n}^{n/2}=\omega_2=-1\)
证明：\(\omega_{dn}^{dk}=(e^{2\pi i/dn})^{dk}=(e^{2\pi i/n})^k=\omega_n^k\)
引理2（折半引理）：如果\(n>0\)为偶数，那么\(n\)个\(n\)次单位复数根的平方的集合就是\(n/2\)个\(n/2\)次单位复数根的集合。
证明：根据消去引理，对任意非负整数\(k\)，有\((\omega_n^k)^2=\omega_{n/2}^k\)。注意，如果对于所有\(n\)词单位复数根进行平方，那么获得每个\(n/2\)次单位根正好\(2\)次，因为\((\omega_n^{k+n/2})^2=\omega_n^{2k+n}=\omega_n^{2k}\omega_n^n=\omega_n^{2k}=(\omega_n^k)^2\)
引理3（求和引理）：对任意整数\(n \geq 1\)和不能被\(n\)整除的非负整数\(k\)，有\(\sum_{j=0}^{n-1} (\omega_n^k)^j=0\)。
证明：
\(k\)不能被\(n\)整除保证分母不为\(0\)。
记
则\(y=(y_0, y_1, ..., y_{n-1})\)就是系数向量\(a=(a_0, a_1, ..., a_{n-1})\)的离散傅里叶变换(DFT)，记为\(y=DFT_n(a)\).
FFT
(以下\(n\)默认为\(2\)的幂)
通过FFT，利用复数单位根的特殊性质，可以在\(O(n log n)\)的时间内算出\(DFT_n(a)\).FFT利用分治策略，采用\(A(x)\)中偶数下标的系数与奇数下标的系数，分别定义两个新的次数界为\(n/2\)的多项式\(A^{[0]}(x), A^{[1]}(x)\):
则\(A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)\)
因此\(A(x)\)在\(\omega_n^0, \omega_n^1, ..., \omega_n^{n-1}\)的值转换为：
- 求次数界为\(n/2\)的多项式\(A^{[0]}(x)\)和\(A^{[1]}(x)\)在点\((\omega_n^0)^2, (\omega_n^1)^2, ..., (\omega_n^{n-1})^2\)的值
- 根据\(A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)\)进行整合
根据折半引理，1中对点求值并不是\(n\)个不同的值，而是\(n/2\)个\(n/2\)次单位复数根。即把一个\(n\)个元素的\(DFT_n(a)\)计算划分为两个规模为\(n/2\)个元素的\(DFT_{n/2}\)计算。
当\(n=1\)时，\(y_0=a_0\omega_1^0=a_0\)
否则设\(y_k^{[0]}=A^{[0]}(\omega_{n/2}^k), y_k^{[1]}=A^{[1]}(\omega_{n/2}^k)\)，根据消去引理，有\(y_k^{[0]}=A^{[0]}(\omega_{n}^{2k}), y_k^{[1]}=A^{[1]}(\omega_{n}^{2k})\)
则
逆运算
用矩阵表示\(y=V_na\)，则\(V_n\)在\((k, j)\)处元素为\(\omega_n^{kj}\)。求出\(V_n^{-1}\)，则可以求出\(a=DFT_n^{-1}(y)\)。
定理：\(V_n^{-1}\)的\((j, k)\)处元素为\(\omega_n^{-kj}/n\)。
证明：考虑\(V_n^{-1}V_n\)中(j, j')处的元素：
若\(j=j'\)，则此和为\(1\)；否则根据求和引理，此和为\(0\)，所以\([V_n^{-1}V_n]\)为单位矩阵。\(a_j=\frac{1}{n}\sum_{k=0}^{n-1}y_k\omega_n^{-kj}\)。
对比DFT：把\(a\)与\(y\)互换，用\(\omega_n^{-1}\)替换\(\omega_n\)，并将结果除以\(n\)。因此逆运算也可以在\(O(nlogn)\)内完成。
注意：
- 不够\(2\)的幂时要补零
- 做完FFT后顺序是乱的，应该按照编号的二进制的翻转对应的十进制进行排序。例如：\(n=8\)时
0 000 000 0
1 001 100 4
2 010 010 2
3 011 110 6
4 100 001 1
5 101 101 5
6 110 011 3
7 111 111 7
右边为对应编号
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <deque>
#include <queue>
#include <vector>
#include <map>
#include <complex>
using namespace std;
typedef complex<double> E;
const int maxn=int(1e6)+100;
const double PI=acos(-1);
int n, m;
int rev[maxn];
E a[maxn], b[maxn];
void init()
{
scanf("%d%d", &n, &m);
for (int i=0, tmp; i<=n; ++i)
{
scanf("%d", &tmp);
a[i]=E(tmp, 0);
}
for (int i=0, tmp; i<=m; ++i)
{
scanf("%d", &tmp);
b[i]=E(tmp, 0);
}
m=n+m;
for (n=1; n<=m; n<<=1);
}
void FFT_init()
{
for (int i=0; i<n; ++i)
for (int j=0; 1<<j<n; ++j)
rev[i]=(rev[i]<<1) | (i>>j & 1);
}
void FFT(E *a, int type)
{
for (int i=0; i<=n; ++i)
if (i<rev[i]) swap(a[i], a[rev[i]]);
for (int i=2; i<=n; i<<=1)
for (int j=0; j<n; j+=i)
{
E w(cos(2*PI/i), sin(type*2*PI/i)), wn(1, 0);
for (int k=0; k<i>>1; ++k, wn*=w)
{
E tmp=a[j+k];
a[j+k]=a[j+k]+wn*a[j+k+(i>>1)];
a[j+k+(i>>1)]=tmp-wn*a[j+k+(i>>1)];
}
}
}
void solve()
{
FFT(a, 1);
FFT(b, 1);
for (int i=0; i<n; ++i) a[i]=a[i]*b[i];
FFT(a, -1);
for (int i=0; i<=m; ++i)
printf("%d ", int(a[i].real()/n+0.5));
}
int main()
{
init();
FFT_init();
solve();
return 0;
}
