FFT—快速傅里叶变换学习小记

我们要干什么

FFT能够快速地解决多项式乘法。额···好像也只能干这个事。
一个特殊的例子：高精度乘法。

问题简述

我们记一个多项式$A(x)$的次数界为n，则
$A(x)=\sum_{i=0}^{n-1}a_i*x^i$，其中a为系数，x为变量。
注意最高次系数为n-1，不是n。实际上次数界就是有多少项系数。
两个多项式相乘，我们一般记为$C(x)=A(x)B(x)$，C的次数界显然是A,B的和-1。一般来说，乘起来要$O(n^2)$

FFT的中心思想

转换思路

普通的相乘方法：$c_j=\sum_{i=0}^{j}a_i*b_{j-i}$
提出概念： 点值，插值。

点值

一个次数界为n的多项式A的点值表达为n个点值对所组成的集合。形如
$\{(x_0,y_0),(x_1,y_1)...(x_{n-1},y_{n-1})\}$
其中$y_i=A(x_i)$，x各不相同。
比如多项式$A(x)=x^2-3x+2$的点值表达可以为，{(-1,6),(0,2),(5,12)}。

插值

插值运算为点值运算的逆运算，如果给定一个n个点值对的点值表达，我们可以确定唯一一个次数界为n的多项式。就比如通过上面的三个点，可以求出A(x)的各项系数。

综合两者进行多项式乘法

确定一组x，可以得到A,B的点值表达。
$A:\{(x_0,y_0),(x_1,y_1)...(x_{n-1},y_{n-1})\}$
$B:\{(x_0,y'_0),(x_1,y'_1)...(x_{n-1},y'_{n-1})\}$
那么可以知道C的点值表达为
$C:\{(x_0,y_0*y'_0),(x_1,y_1*y'_1)...(x_{n-1},y_{n-1}*y'_{n-1})\}$
$不太理解的话，可以把y_0和y'_0的表达式写出来相乘，发现y_0*y'_0的表达式和前者确实相同$
通过插值运算，就可以得到多项式C的系数了。
注意：A,B次数界均为n的话，C的次数界为2n-1(为了方便直接弄成2n)，所以x集合大小至少为2n，不然插值失败。

流程

对于次数界均为n的两个多项式A,B
1.点值运算，构造长度为2n的点值表达；
2.每个点的y相乘，算出C的点值表达；
3.插值运算，通过点值表达求出系数。
然而暴力来做时间复杂度仍为$O(n^2)$，没有区别。
其实我们可以通过选特殊的x集合，来优化时间，我们选的数叫做n次单位复数根。

N次单位复数根及其性质

n次单位复数根是满足$w^n=1$的复数$w$。那么我们知道给定了n，n次单位复数根恰好有n个。对于根k=0,1…n-1，这些根是$e^{2\pi ik/n}(i为虚数单位)$，为了解释，我们需要用到复数的指数形式的定义。
$e^{iu}=cos(u)+isin(u)$。
这一条是欧拉定理(之一)，实际上，我们发现$e^{iu}$的泰勒展开和$cos(u)$的泰勒展开加上$isin(u)$的泰勒展开是一样的。所以两边相等。
我们来看个图，直观地发现规律。

可以发现，他们均匀分布在复平面(以实数域为x轴，虚数域i为y轴的平面)的单位圆上。跟学三角函数的时候那些东西是一样的。

性质1，主n次单位根

我们称$w_n=w_n^1$为主n次单位根，注意到每个n次单位复数根都是由旋转得到的，每次旋转有一个固定的角度，可见$w_n^i=w_n^{i-1}*w_n$。
n次单位复数根可视为公比为$w_n$的等比数列。

性质2，群的性质

由于单位根们可表达为一对三角函数，那么他也是具有周期的
$w_n^0=w_n^n=1,w_n^k=w_n^{k\%n}$
而由性质1得，一个$w_n^k$实际上可以看成$(w_n)^k$，那么$w_n^k*w_n^j=w_n^{j+k}$
又可以得到$(w_n^k)^2=w_n^{2k}$

两个引理

消去引理

$w_{dn}^{dk}=w_n^k$，可以从图象角度理解，他们所对应的点是相同的。
推论：
$w_n^{n/2}=w_2=-1$
$w_n^{2k}=w_{n/2}^k$

折半引理

若n为>0的偶数，那么n次单位复数根的平方的集合，等于n/2次单位复数根的集合。
这个也好理解，可以知道$(w_n^k)^2=(w_n^{k+n/2})^2$，这个由上面的一些性质可以推出来。图象上理解：n次单位复数根的平方，就对应着三角函数里的角度乘了2。

求和引理

对于任何整数n>=1和不能被n整除的非负整数k（负数可以通过群的性质变为正数），有
$\sum_{j=0}^{n-1}(w^k_n)^j=0$
等比数列求和，化出来发现末项和首项都是1，所以是0。k为n的时候分母为0，所以不成立。

FFT及其关键算法

DFT——离散傅里叶变换

我们的目标：计算次数界为n的多项式$A(x)=\sum_j a_j*x^j$在n次单位复数根上的值。
定义结果y
$y_k=A(w_n^k)=\sum_{j=0}^{n-1}a_j*w_n^{jk}$
y即为a的离散傅里叶变换。
可记为$y=DFT_n(a)$

FFT——快速傅里叶变换

为了充分利用n次单位复数根，我们令n为2的幂，这样可以在$O(nlog_2n)$时间内求出$DFT_n(a)$。为了这样做，我们用分治的策略。

分治策略

如何求出单个数x的函数值A(x)？我们定义两个多项式
$A_0(x)=a_0+a_2x+a_4x^2···a_{n-2}x^{n/2}$
$A_1(x)=a_1+a_1x+a_3x^2···a_{n-1}x^{n/2}$
即一个全是偶数编号的系数，一个全是奇数的。我们发现这两个东西的次数是n/2，变小了一半。
可得$A(x)=A_0(x^2)+xA_1(x^2)$，那么再使用点值相乘。

阶段性总结

原问题：求$A(x)$在n次单位复数根上的函数值。
转化成：求$A_0(x)和A_1(x)$在n/2次单位复数根上的值。
可以发现不同的自变量少了一半，而每个自变量出现2次。
于是，我们递归地对n/2的多项式$A_0(x)$与$A_1(x)$在n/2个n/2次单位复数根进行求值。

伪代码

看到最后的处理，我们对于$(w_n^k)^2$相同的一起处理，而平方相同的两个n次单位复数根恰好互为相反数。
运行时间显然是$O(nlogn)$的。
这样我们解决了点值运算。可以处理出C的点值表达了！

插值运算

这个其实跟点值运算差不多。
我们如果把点值运算写成矩阵方程的形式，可以得到表达式$y=V_na$。

只是把$\sum$换个形式写而已。
那么可以知道$a=y*V_n^{-1}$。也就是说，只要得到逆矩阵，再像上面做类似的点值运算，就能得到系数了！

逆矩阵

定理，$V_n^{-1}[j,k]=w_n^{-jk}/n$。
证明很麻烦，考虑证明$V_n*V_n^{-1}=单位矩阵$，即(i,i)处为1，其他为0的矩阵。
证明：
$[V_n*V_n^{-1}]_{j,j'}=\sum_{k=0}^{n-1}w_n^{k*(j'-j)}/n$，如果j’=j，那后面为1，否则根据求和引理，为0。
那么，只需要把原来fft里面的a和y互换，$w_n变成w_n^{-1}$，再进行fft，结果再除n，就能得到系数了！我们称这个为逆DFT。

实现的优化

我们不可能像伪代码一样赋值数组嘛···考虑给数组重新分配位置，使得最后算出来又是原来的顺序。
观察这幅图

最下面一层为我们数组应该保存的顺序。可以发现，他们的位置，刚好就是原下标在二进制下的翻转(开头为0，结尾为次数界n在二进制下的最高位的翻转)。
然后对于上面几层，我们每个组里面有多个值。

发现如果按照刚刚那个顺序存，每一层转移到下一层的时候，一个$A(x_0)$要提取的$A_0(x_0)和A_1(x_0)$的间隔都是一样的。
设一行分为cnt组，即相同的$x_0$出现cnt次，共有m个不同的$x_0$，可以发现转移间隔为m/2。而平方相等的n次单位复数根的间隔也是m/2，这就很优美了。
这个大概理解一下就好，反正代码很简洁明了。

代码

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<bitset>
#include<cmath>
#define fo(i,j,k) for(i=j;i<=k;i++)
#define fd(i,j,k) for(i=j;i>=k;i--)
typedef long long ll;
typedef double db;
using namespace std;
const int N=400005;
const db pi=acos(-1),eps=1e-5;
struct vec
{
    db a,b;
    vec (db ax=0,db bx=0) {a=ax,b=bx;}
};
vec operator +(vec a,vec b)
{
    return vec(a.a+b.a,a.b+b.b);
}
vec operator -(vec a,vec b)
{
    return vec(a.a-b.a,a.b-b.b);
}
vec operator *(vec a,vec b){return vec(a.a*b.a-a.b*b.b,a.a*b.b+a.b*b.a);}
int Log,m,len,i,n,ans[N];
vec t[N],a[N],b[N];


void DFT(vec *a,int n,int sig)
{
    int i,j,k,half,siz;
    vec u,v,w;
    fo(i,0,n-1)
    {
        int p=0;
        for(int pos=0,tmp=i;pos<Log;pos++,tmp/=2) p=(p<<1)+(tmp&1);
        t[p]=a[i];
    }
    for(int m=2; m<=n; m*=2)
    { 
        int half = m/2; 
        for(int i=0; i<half; i++)
        { 
            vec w( cos(i*sig*pi/half), sin(i*sig*pi/half));
            for(int j=i; j<n; j+=m)
            { 
                k=j+half;
                u=t[j];v=w*t[k];
                t[j]=u+v;
                t[k]=u-v;
            }
        }
    }

    fo(i,0,n-1) a[i]=t[i];
}

void FFT(vec *a,vec *b)
{
    DFT(a,len,1);
    DFT(b,len,1);
    fo(i,0,len-1) a[i]=a[i]*b[i];
    DFT(a,len,-1);
    fo(i,0,len-1) a[i].a/=(db)len;
}
int max2(int x)
{
    int ret;
    for(ret=1;ret<x;ret<<=1);
    return ret<<1;
}
int ci(int x) {return (x==1)?0:ci(x/2)+1;}
int main()
{
    freopen("taxi.in","r",stdin);
    scanf("%d %d",&n,&m);
    n++;m++;
    len=max2(max(n,m));
    Log=ci(len);
    fo(i,0,n-1) scanf("%lf",&a[i].a);
    fo(i,0,m-1) scanf("%lf",&b[i].a);
    FFT(a,b);
    fo(i,0,n+m-1) ans[i]=round(a[i].a+eps);
    fo(i,0,n+m-2) printf("%d ",ans[i]);
}
//注意，由于是实数运算，所以给出系数在10^5以上时很容易出现误差

总结

中心思想：点值运算与插值运算
优化时间复杂度：选取n次单位复数根
难点：n次单位复数根的各种性质
重点：FFT的分治思想，DFT与逆DFT
要背的地方：FFT实现优化