背板子.jpg
快速傅里叶变换(FFT) 学习笔记
特别鸣谢
以及
什么是FFT
什么是多项式乘法
现在加入我们有两个多项式：
$f(x)=a1x^2+b1x+c1$
$g(x)=a2x^2+b2x+c2$
那么如果暴力把这两个多项式相乘，就是lena*lenb的复杂度，也就是O(n^2)的
那FFT呢
我们可以先理解为，FFT是一种可以把O(n^2)的多项式乘法转化为O(n log n)的算法
前置知识
系数表示法和点值表示法
系数表示法就是我们平时用的表示法，用n+1个系数表示一个n次函数
$f(x)=a0+a1x+a2x^2+..+anx^n ⇔ f(x) = (a0,a1,a2,..,an) $
点值表示法就是用n+1个点来表示一个n次函数
然后通过高斯消元解方程组可以得到原函数
DFT 和 IDFT
把一个多项式转化为离散的点的方法就叫做DFT(离散傅里叶变换)
把一堆离散的点还原回多项式的方法就叫做IDFT(离散傅里叶反变换)
FFT(快速傅里叶变换)就是一种用来计算DFT和IDFT的快速算法
虚数
相对于“实数”的概念
如果说$\sqrt{1}$是一个实数单位的话，那么我们定义$\sqrt{-1}$是一个虚数单位
复数
首先咱们定义一个平面坐标体系，叫做复平面
这个复平面以实数部为x轴，以虚数部作为y轴，平面上任意个点称作极坐标
设某个极坐标为(a,b)，与原点的连线和x轴的夹角为θ，长度为r
那么这个点除了(a,b)还可以用a+bi或(r,θ)表示，这就是一个复数
我们可以认为，两个复数相乘有特殊的几何意义，即
$(r1,θ1)(r2,θ2)=(r1r2,θ1+θ2)$
单位复根
如果我们有两个用点值表示的多项式
$f(x)={(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),..,(x_n,f(x_n))}$
$g(x)={(x_0,g(x_0)),(x_1,g(x_1)),(x_2,g(x_2)),…,(x_n,g(x_n))}$
那么这两个多项式相乘就是
$f(x)\times g(x)={(x_0,f(x_0)g(x_0)),(x_1,f(x_1)g(x_1)),(x_0,f(x_2)g(x_2)),..,(x_n,f(x_n)g(x_n))}$
染鹅
这样复杂度贼TM高啊
但是我们可以通过找到一个x值，使得我们不用运算n次方操作。
那么实数就有1和-1，虚数就有i和-i
可这只有4个啊
我们考虑复数，在复平面上，以原点为圆心，画一个半径为1的圆，那么圆上的点相乘时，r永远等于1，只是 θ的大小在发生改变，而这个圆上有无穷个点能满足它的某次方等于1。
我们把这种数叫做“复根”，用ω表示，如果$ω^k=1$，那么我们称“ω为1的k次复根”计做“$ω_k^n$”
其中n就是单位复根的n次方咯——那么，$ω_k^1$就是单位复根本身了
当然你可能会问这玩意怎么用程序写出来——这就要看代码了
FFT的实现
DFT
运用了分治的思想来求解，分治多项式的“奇数次项”和“偶数次项”。
$f(x)=a_0+a_1x+a_2x^2+…+a_7x^7$
$f(x)=(a_0+a_2x^2+…+a_6x^6)+x*(a_1+a_3x^2+…+a_7x^6)$
那么令
$g(x)=(a_0+a_2x+…+a_6x^3)$
$h(x)=(a_1+a_3x+…+a_7x^3)$
则
$f(x)=g(x^2)+x*h(x^2)$
这个二分也是有局限性的，就是它的长度不一定是2的整数倍，我们就先把它补成2的整数倍的长度，当然，最高项的系数为0。
我们还可以对这个分治的过程继续优化，因为如果多项式足够长，递归的过程要消耗大量的内存，我们就可以先在原数组内“拆分”，然后用一种类似倍增的方法计算多项式的值。
我们先看原先的分治过程是怎样的
$a_0,a_1,a_2,a_3,…,a_7$
$(a_0,a_2,a_4,a_6),(a_1,a_3,a_5,a_7)$
$(a_0,a_4),(a_2,a_6),(a_1,a_5),(a_3,a_7)$
$a_0,a_4,a_2,a_6,a_1,a_5,a_3,a_7$
我们把原数组和分治之后的对比
$a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7$
$a_0,a_4,a_2,a_6,a_1,a_5,a_3,a_7$
这样还看不出什么来，可如果我们把下标转换成二进制呢？
$a_{000},a_{001},a_{010},a_{011},a_{100},a_{101},a_{110},a_{111}$
$a_{000},a_{100},a_{010},a_{110},a_{001},a_{101},a_{011},a_{111}$
仔细一看，这不就是把二进制位从左到右翻转过来吗！
这里就不证明了……
IDFT
这个东西……意会一下吧
笔者实在是难以详尽地理解或者描述IDFT的流程 推荐Blog
列一下公式
我们已知 $e^{iπ}=-1$
则 $e^{2iπ}=1$
$w_k=e^{2iπ/k}$
$1/w_k=e^{2i*(-π)/k}$
然后结合底下的代码理解一下
代码实现
#include<iostream>
#include<cstdio>
#include<cmath>
#define N (4000000+100)
#define pi acos(-1)
using namespace std;
int n,m,fn,len,r[N];
struct cop {
double x,y;
cop (double xx=0,double yy=0) {
x=xx; y=yy;
}
}a[N],b[N];
cop operator + (cop a,cop b) {return cop(a.x+b.x,a.y+b.y);}
cop operator - (cop a,cop b) {return cop(a.x-b.x,a.y-b.y);}
cop operator * (cop a,cop b) {return cop(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
cop operator / (cop a,double b) {return cop(a.x/b,a.y/b);}
void FFT(int n,cop *a,int opt) {
for (int i=0; i<n; i++) if (i<r[i]) swap(a[i],a[r[i]]);
for (int k=1; k<n; k<<=1) { //当前分治段的长度
cop wn=cop(cos(pi/k),opt*sin(pi/k));
for (int i=0; i<n; i+=(k<<1)) { //当前分治段的起点
cop w=cop(1,0);
for (int j=0; j<k; j++,w=w*wn) { //枚举段内元素
cop x=a[i+j],y=w*a[i+j+k];
a[i+j]=x+y; a[i+j+k]=x-y;
}
}
}
if (opt==-1) for (int i=0; i<n; i++) a[i]=a[i]/n;
}
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);
fn=1;
while (fn<=n+m) fn<<=1,len++;
for (int i=0; i<fn; i++)
r[i]=(r[i>>1]>>1) | ((i&1)<<(len-1));
FFT(fn,a,1); FFT(fn,b,1);
for (int i=0; i<=fn; i++) a[i]=a[i]*b[i];
FFT(fn,a,-1);
for (int i=0; i<=n+m; i++)
printf("%d ",(int)(a[i].x+0.5));
}
注意事项
假如说我有一个点值表达的多项式，我能不能直接用IDFT把它转化成系数表达式？
不能！
因为IDFT是在DFT的基础上进行的操作，也就是说IDFT中的“点值”是虚数形式的特殊点值，换成正常的点值表达式来操作是不会正确的。
