快速傅里叶变换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^{2pi ik/n}, k=0, 1, 2, ..., n-1, omega_n=e^{2pi 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^{2pi i/dn})^{dk}=(e^{2pi 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_0omega_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_komega_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;
}