【算法】FFT 食用指北 (1) —— by ScapeGoatTree

引言

占了好长时间的位置，思来想去终于打算写这篇文章了 qwq

关于 FFT，印象最深的是 2017 在 NOIP 的前夜，和一群大佬们聚众狼人杀的时候，我校某省队爷放出的话：FFT 这个东西，在我们城市，你想学，有两种方法，第一，你来问我们三个进过省队的；第二，你自学。

于是我选择了自学。

不知道看了多少 dalao 们多少不同思路的博客，也不知道看过多少种不同花式写法的 FFT 代码，总算自己囫囵吞枣背得一个板子，写得几个题目，便也大言不惭地来胡乱写写指北了。不为别的，我不希望当其他人想学的时候无人可问

简介

FFT，快速傅里叶变换，是一种在 $O(nlog_2n)$的复杂度内解决 DFT（离散傅里叶变换）常规算法 $O(n^2)$复杂度解决的问题。DFT 本身作用我们并不需要了解，我们需要知道的是 FFT 可以做什么。

FFT 可以用来求多项式卷积。

什么叫多项式卷积？

对于一个多项式 $F(x)=\sum\limits _ {i=0}^n a _ ix^i$和 $G(x)=\sum\limits _ {i=0}^m b _ ix^i$, 我们定义卷积 $C(x)=F(x) * G(x)=\sum\limits _ {i=0}^{n+m} \sum\limits _ {j=0}^i a _ j * b _ {i-j}x^i$

也就是说，$C(x)=a _ 0b _ 0+(a _ 1b _ 0+a _ 0b _ 1)x+(a _ 2b _ 0+a _ 1b _ 1+a _ 0b _ 2)x^2….$

观察一下，就会发现其实就是两个多项式各项相乘得出来的结果。
给定两个多项式的情况下，如果我们想要计算其卷积的结果，那么我们通常会采取一种方法：先取值，后插值。取出 n 个点，计算他们在多项式中的值，分别相乘，这个结果一定是新多项式在这个点的值。时间复杂度常常是 $O(n^2).$

引例

先看一个例题：FFT 模板题
这明确要求我们求卷积。n 是 $10^6$范围，我们用普通的想法是解决不了的。这个时候，我们就要开始介绍 FFT 了。

FFT 算法的宗旨是，通过引入复数单位根的方法，在瓶颈步骤加以速度优化。那么现在我们必须找到瓶颈步骤在哪里。

取值要进行计算。取点和计算的时间复杂度共计是 $O(n^2)$的。牛顿插值，泰勒展开和拉格朗日插值法都是 $O(n^2)$的时间复杂度。这个复杂度我们接受不了。

两方面都是大问题，那我们本着逐个击破的原则先解决取点的问题。

剖析

我一定要在实数域内取点子吗？

正常情况下我们处理实数多项式不可能牵扯到实数域外的世界。但是这次，复数单位根帮了我们一个大忙。

实数域内的数字很高傲，一是一二是二，取点 1 可以很快地得出结果，但是 2，3 就不行。他们各自直接互不相同。但是单位根不太一样，他们之间有着千丝万缕的联系。

在复平面上，以原点为圆心，1 为半径作圆，所得的圆叫单位圆。以圆点为起点，圆的 n 等分点为终点，做 n 个向量，设幅角为正且最小的向量对应的复数为 $\omega_n$，称为 n 次单位根。

很容易就可以看出来，$\omega_n^0=\omega_n^n=1$, 但是其他单位根我们不好解决。怎么办？有一座桥梁，沟通了虚数和实数。

• 欧拉公式 (应用于单位根)：$\omega _ n^k=cos(k * \frac{2\pi}{n})+isin(k * \frac{2\pi}{n})$

在欧拉公式的帮助之下，我们克服了虚数求值的问题。从几何关系来看，我们还可以看出如下两个性质：
– $\omega_n^k=\omega_{2n}^{2k}$ 这两个矢量是重合的。
– $\omega_n^k=-\omega_n^{k+\frac{n}{2}}$ 这等价于绕了 180°，关于原点对称自然是负号

接下来我们将 FFT 的公式做一个详细地推导。至少，我们要从原理上理解这个算法。

求值

考虑有一个多项式 $F(x)=a_0+a_1x+a_2x^2+…+a_nx^n$
现在我们将他的偶次项分离，得到两个多项式。为了方便，我们假定 n 是奇数：
$F _ 1(x)=a _ 1x+a _ 3x^3+…+a _ {n}x^{n}$和 $F _ 2(x)=a _ 0+a _ 2x^2+…+a _ {n-1}x^{n-1}$

令 $G _ 1(x)=a _ 1+a _ 3x+…+a _ nx^{\frac{n}{2}-1}$和 $G _ 2(x)=a _ 0+a _ 2x+…+a _ {n-1}x^{\frac{n}{2}-1}$

我们可以发现 $F(x)=G _ 2(x^2)+xG _ 1(x^2)$
为了求出多项式 $F(x)$的系数，我们需要将所有的单位根带进去分别求值。带入 $\omega _ n^k$和 $\omega _ n^{k+\frac{n}{2}}$时，我们发现
$F _ k=G _ 2(\omega _ n^{2k})+\omega _ n^kG _ 1(\omega _ n^{2k})$, 并且此时 $F _ {k+\frac{n}{2}}=G _ 2(\omega _ n^{2k})-\omega _ n^kG _ 1(\omega _ n^{2k})$
如果我们分别求出了 $G _ 1,G _ 2$在 $\omega _ {2n}^k$时的值，我们就可以 $O(1)$地求出两个点代表的值。于是我们可以缩小一半的规模。

缩小一半就结束了？不，没有结束。注意到 $G _ 1,G _ 2$是从原来多项式中分离出来的，他们本身也是未知的。好办，我们对于两个小问题分别使用 FFT 再进行一次求值。这样做，我们每次可以将问题的规模缩小一半。

运用主定理，我们可以得出 FFT 算法求值的时间复杂度变为了 $O(nlog _ 2n)$
但是光求值优化不管用，现在我们必须插值回去来确定这个多项式。

插值

整个求值过程中，我们都在将一个特定的值——单位根不断地带入他的各个幂中，这令我们想到了范德蒙德矩阵：
$V(x _ 0,x _ 1,\cdots ,x _ {n-1})=\prod\limits _ {n > i > j \geq 0}(x _ {i}-x _ {j})$
这是他的通项公式，这家伙长成这样：
$V=\begin{pmatrix} {1}&{1}&{\cdots}&{1}\\{x _ {0}}&{x _ {1}}&{\cdots}&{x _ {n-1}}\\{x _ {0}^2}&{x _ {1}^2}&{\cdots}&{x _ {n-1}^2}\\{\vdots}&{\vdots}&{}&{\vdots}\\{x _ {0}^{n-1}}&{x _ {1}^{n-1}}&{\cdots}&{x _ {n-1}^{n-1}}\\\end{pmatrix}$

如果我们构造一个向量 $\begin{pmatrix} a _ 0 \\a _ 1 \\a _ 2 \\ \cdots\\a _ {n-1} \\ \end{pmatrix}$就会发现 $\begin{pmatrix} a _ 0 \\a _ 1 \\a _ 2 \\ \cdots\\a _ {n-1} \\ \end{pmatrix}*\begin{pmatrix} {1}&{1}&{\cdots}&{1}\\{x _ {0}}&{x _ {1}}&{\cdots}&{x _ {n-1}}\\{x _ {0}^2}&{x _ {1}^2}&{\cdots}&{x _ {n-1}^2}\\{\vdots}&{\vdots}&{}&{\vdots}\\{x _ {0}^{n-1}}&{x _ {1}^{n-1}}&{\cdots}&{x _ {n-1}^{n-1}}\\\end{pmatrix}=\begin{pmatrix} y _ 0 \\y _ 1 \\y _ 2 \\ \cdots\\y _ {n-1} \\ \end{pmatrix}$,

其中 y 都是刚才带入的各个点值 x 求出来的结果。这就表明，一个表示多项式系数的向量乘上范德蒙德矩阵就会相当于做了一次多项式多点求值。那么我们只需要将手上的 y 向量呈上范德蒙德矩阵的逆矩阵就可以了。

逆矩阵怎么构造呢？这部分内容建议去看线性代数，他讲的肯定比我清楚。当然不看也可以，因为对于 FFT 我们要掌握他的原理，数学推导只需要注意几个重要的地方就行。这里直接给出结果：
$V(w _ 0,w _ 1,\cdots,w _ {n-1})^{-1}={1 \over n}V(w _ 0,w _ {-1},\cdots,w _ {-n+1})$
于是我们就可以将单位根全部取负，将这些新的点带入做一次 FFT 求值就可以了 qwq。得到的结果乘上（1/n）就得到了系数！

实现

注意到我们每次会将问题一分为二，为了防止繁杂的边界讨论，我们统一将最高位补足到 2 的整数次幂。
根据之前的实现原理，我们可以比较容易的写出如下的代码：

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(cos(2.0 * Pi / limit) , type * sin(2.0 * Pi / limit)), w = complex(1, 0);
    //Wn 为单位根，w 表示单位根幂
    for (int i = 0; i < (limit >> 1); i++, w = w * Wn) //这里的 w 相当于公式中的 k
        a[i] = a1[i] + w * a2[i],
               a[i + (limit >> 1)] = a1[i] - w * a2[i]; //利用单位根的性质，O(1) 得到另一部分
               //这里不妨说一下，我的代码是参照 attack 的模板学习的
}

但是递归操作涉及到大量的栈操作，于是 Fast Fourier Transformation 变成了 Fast Fast TLE。

那么，通过一种巧妙的数学关系——二进制倒序重排，我们可以发现我们求和的顺序是符合一个重排之后的 “递增序列” 的。这就免去了我们在计算时先按照奇偶系数排列的事情了。取而代之的是 rader 操作，我们能够快速地找出每个位置所对应的 $O(1)$可求的对象位置。现在，我们可以写出一个更加漂亮的迭代实现的 FFT 算法：

inline void FFT(comp A[],int mode)
{
    for(register int i=0;i<lim;++i)
        if(i<rader[i]) swap(A[i],A[rader[i]]);
    for(register int mid=1;mid<lim;mid<<=1)
    {
        comp Wn;
        Wn.x=cos(Pi/mid),Wn.y=mode*sin(Pi/mid);
        for(register int R=mid<<1,j=0;j<lim;j+=R)
        {
            comp w(1,0);
            for(register int k=0;k<mid;++k,w=w*Wn)
            {
                comp x=A[j+k],y=w*A[j+mid+k];
                A[j+k]=x+y;
                A[j+mid+k]=x-y;
            }
        }
    }
}

其中 rader 数组用于 rader 操作，我们可以按照这样的方法提前处理好：

for(register int i=0;i<lim;++i)
        rader[i]=(rader[i>>1]>>1)|((i&1)<<(l-1));

复数的运算可以使用 STL 中的 complex 模板，如果您是 STL 重度恐惧症患者，您可以按照如下方式定义自己的模板：

struct comp
{
    double x,y;
    comp (double xx=0,double yy=0) {x=xx,y=yy;}
};
comp operator + (comp a,comp b) {return comp(a.x+b.x,a.y+b.y);}
comp operator - (comp a,comp b) {return comp(a.x-b.x,a.y-b.y);}
comp operator * (comp a,comp b) {return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

贴出完整版的代码。请不要直接使用这个代码进行提交，因为这会导致你失去熟练他的机会。

#include<bits/stdc++.h>
#define maxn 10000010
using namespace std;
inline int read()
{
    int n=0,k=1;
    char c=getchar();
    while(c>'9'||c<'0') {if(c=='-') k=-1; c=getchar();}
    while(c<='9'&&c>='0') n=(n<<1)+(n<<3)+(c^48),c=getchar();
    return n*k;
}
const double Pi=acos(-1.0);
struct comp
{
    double x,y;
    comp (double xx=0,double yy=0) {x=xx,y=yy;}
}acid[maxn],base[maxn];
comp operator + (comp a,comp b) {return comp(a.x+b.x,a.y+b.y);}
comp operator - (comp a,comp b) {return comp(a.x-b.x,a.y-b.y);}
comp operator * (comp a,comp b) {return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
int l,rader[maxn];
int lim=1;
inline void FFT(comp A[],int mode)
{
    for(register int i=0;i<lim;++i)
        if(i<rader[i]) swap(A[i],A[rader[i]]);
    for(register int mid=1;mid<lim;mid<<=1)
    {
        comp Wn;
        Wn.x=cos(Pi/mid),wn.y=mode*sin(Pi/mid);
        for(register int R=mid<<1,j=0;j<lim;j+=R)
        {
            comp w(1,0);
            for(register int k=0;k<mid;++k,w=w*Wn)
            {
                comp x=A[j+k],y=w*A[j+mid+k];
                A[j+k]=x+y;
                A[j+mid+k]=x-y;
            }
        }
    }
}
int n,m;
int main()
{
    n=read(),m=read();
    for(register int i=0;i<=n;++i) acid[i].x=read();
    for(register int i=0;i<=m;++i) base[i].x=read();
    while(lim<=n+m) lim<<=1,++l;
    for(register int i=0;i<lim;++i)
        rader[i]=(rader[i>>1]>>1)|((i&1)<<(l-1));
    FFT(acid,1);
    FFT(base,1);
    for(register int i=0;i<=lim;++i) acid[i]=acid[i]*base[i];
    FFT(acid,-1);
    for(register int i=0;i<=n+m;++i)
        printf("%d ",(int)(acid[i].x/lim+0.5));
        //0.5 是为了避免精度误差
    printf("\n");
    return 0;
}

从而，我们也圆满地完成了 FFT 算法的学习。建议将这个算法背出来。目前我所能达到的水平是 100 分钟完成默写+调试+ac。至少不能高于 20 分钟，否则会耽误你大量的考场时间。

以上就是这篇文章的全部内容，更多关于多项式算法的知识我将会在之后的文章中给出。