FFT详细教程

阅读原文时间：2021年04月22日阅读：1

多项式简介

什么是多项式

众所周知，多项式可以表示为：
\[f(x)=\sum ^{n} _{i=0} a_i x^i\]
其中，\(a_i\)表示为多项式的系数。

还有一种表示方法是用一个\(n+1\)维的向量表示这个多项式，我们称为该多项式的系数表达。如下：
\[\vec{a}=(a_0, a_1, \cdots, a_n)\]
系数表达在字面上就能知道这个表达的意思了。

而它们的最高非零项的次数称为多项式的度数，用\(\deg\)表示，比如度数为\(n\)的多项式的度数计为\(\deg f(x)=n\)。

多项式乘法

我们定义两个多项式：
\[f(x)=\sum ^{n} _{i=0} a_i x^i\]
\[g(x)=\sum ^{n} _{i=0} b_i x^i\]
则它们的乘法就定义为：
\[h(x)=f(x)\times g(x)=\sum ^{2n} _{i=0} \sum ^{i} _{j=0} a_i b_{i-j} x^i\]
\[\forall \ i>n,\ a_i=b_i=0\]
另外，对于其系数表达做的“乘法”成为卷积，在算法设计中我们一般只讨论卷积，毕竟大多时候只考虑系数。系数的关系为：
\[c_i=\sum^{i}_{j=0}a_jb_{i-j}\]
卷积一般记作：
\[c=a\otimes b\]
那么在后面我们可能将用卷积替代多项式相乘。

像上面那种和式，在算法竞赛中的多项式一类的题目中一般是作为最后推出来的式子的，其特征就是相乘的系数的下标和为定值。

我们发现这个朴素的计算卷积的时间复杂度是\(O(n^2)\)的，那有没有更优的算法使得计算卷积的速度更快呢？

答案是有的，下面我们将由浅入深讲解这个加速算法——快速傅里叶变换算法（FFT）。

当然，在此之前，我们得换一种求多项式的方法。我们可能会绕一个弯子，但是最后的结果会令你大吃一惊。

绕一个弯子

先总体概括一下，这个“弯子”的主要步骤是：先将两个待卷积的多项式分别表示为\(n+m+1\)个不同的点确定的“函数”，然后再将同横坐标的点的纵坐标值相乘，得到确定\(h\)的一些点，再用点值算出\(h\)的系数即可。

Although we might think it is a little complex, we will see why it works in the end. Trust me. It really works.

点值表达

看了总体概括，我们可能还是有点摸不着头脑（真叫人头大），不过别担心。想想我们熟悉的一次函数和二次函数，我们还记得“两点确定一条直线”和“三个不共线的点确定一个抛物线”，运用类比的思维，我们就能理解为什么一个\(\deg f=n\)的多项式可以用\(n+1\)个不同的点来确定了。不过严谨的说法还要运用到范德蒙德行列式来说，这里大家感性地理解就好了。

而我们需要\(n+m+1\)个点是因为多项式\(h\)需要这么多点才能确定。

那么求点值的算法呢？暂时不谈算法优化，我们的求法就是“秦九韶算法”。将\(x\)一个一个提出，然后就可以在复杂度\(O(n+m)\)内可以求出\(1\)个点值，求出我们需要的点值的话，复杂度就是\(O((n+m)^2)\)了。

怪不得是绕了弯子，才第一步就达到了前面我们提到的朴素求法的时间复杂度，要是我走这条路，我多半就会放弃这个法子了（OwO）。

第二步我就不开新的小节了，按照总体概括中提到的，求出\(h\)的点值，时间复杂度仅\(O(n+m)\)。

拉格朗日插值法

有了点值，我们就可以确定多项式了，但是我们似乎不好怎么求\(h\)的系数。其实很简单，按照\(x_i\)对应\(y_i\)我们可以构造出多项式：
\[L(x)=\sum^{n+m+1}_{i=0}\frac{\prod_{i\neq j}(x-x_j)}{\prod_{i\neq j}(x_i-x_j)}y_i\]
虽然看起来十分复杂，但是稍微验算一下我们就知道了。当\(x=x_i\)的时候，我们看到和式中除了一项没有乘以\((x-x_i)\)的（因为式子中要求的\(i\neq j\)），其它的都因该项变为了\(0\)而导致整体被去掉，而剩下的这一项中，分子分母进行约分，最终便剩下了\(y_i\)。

这种插值法的时间复杂度呢？依然是\(O((n+m)^2)\)。

至此，我们也走完了这个“弯子”了。确实不愧叫“弯子”啊，复杂度都多了几倍！这肯定不是我们想要的快速算法。于是我们就需要进行一些优化。

不过我们发现第二步并不需要优化（Even the blind can see it），所以我们便想办法优化第一步和第三步。

这一优化方法令人拍手叫好，让我们一起来看看吧！

快速傅里叶变换算法

算法梗概

该算法优化的方式十分优秀，它并不是对那些求和式进行一大堆复杂的变化，而是对\(x\)的取值进行微妙的操作。在此之前，我们需要引入单位复数根，因为它有一些很好的性质，以达到快速的目的。最后我们再利用矩阵的一些性质重构出我们的结果多项式。

请注意：前方一波操作猛如虎，请站稳扶好。如果有些同学没有学过一些必要的知识，大家可以在文章结尾查看。

单位复数根简介

欧拉公式

微妙操作的起源——欧拉公式。不得不说，欧拉真是数学天才，论文发表狂热者2333。

不扯别的了，我们一起来看看欧拉公式吧。
\[e^{i\theta}=\cos\theta+i\sin\theta\]
至于这是怎么推导出来的，就涉及到一些高数的知识了，我们这里暂且不提，记住就好了。

它的几何意义就是一个处于复平面的圆上的一个弧度为\(\theta\)的点。

进一步的操作

注意到\(\theta=2k\pi\)的情况，我们有\(e^{2k\pi i}=\cos2k\pi+i\sin2k\pi=1\)

考虑\(x^{n}=1\)的解，我们得到\(x=e^{\frac{2k\pi i}{n}}\)，而\(e^{\frac{2\pi i}{n}}\)就是我们要的单位复数根，字面意思对吧！我们一般将其记为\(\omega_{n}\)。取其\(0\)到\(n-1\)的幂次得到的这些值称为\(n\)次单位根。

具象理解：把单位圆劈成\(n\)份就可以了，很形象吧。

单位复数根的性质

我们要进行快速计算的话，没有一些优美有趣（令人窒息）的性质的话，是达不到我们的目的的。而这些性质需要一个特殊的\(n\)的取值：\(2^{k}\)（这和后面涉及到的二分操作有关），这个取值在上面提到的具象理解看来，可以算得上完美分配了（说白了就是强迫症，强迫症万岁）。

消去引理

\[\omega^{dk}_{dn}=\omega^{k}_{n}\]
这个很显然，拿它原来的样子代进去就明白了。

另一个我们需要的就是：
\[\omega^{k+n}_{2n}=-\omega^{k}_{2n}\]
证明的话稍微推一下就好了：
\[\omega^{k+n}_{2n}=e^{\frac{k\pi i}{n}+\pi i}=e^{\frac{k\pi i}{n}}(\cos\pi +i\sin\pi)\]
而\(\cos\pi +i\sin\pi\)求出来就是\(-1\)，再用习惯记法写回来就行了。

折半引理

\(n\)次单位根的平分的集合恰好是\(\frac{n}{2}\)次单位根的集合。即：
\[\omega^{2k}_n=\omega^{k}_\frac{n}{2}\]
这个理解起来也很简单，根据消去引理来推就可以了。

还有重要的一个（和后面的二分计算有关）：
\[(\omega^{k+\frac{n}{2}}_{n})^2=(\omega^{k}_{n})^2\]
证明的话不难，把消去引理的第二个式子的\(n\)都除以\(2\)，然后两边平方就得到了。

虽然这都是根据消去引理推来的，但是记住它不亦乐乎？毕竟我们后面也要用到。

求和引理

\[ \sum^{n-1}_{i=0}(\omega^{k}_{n})^i= \left\{\begin{matrix} 0,n\nmid k \\n,n\mid k \end{matrix}\right. \]
对于\(n\mid k\)的情况，很明显里面全都变成了\(1\)，这个不难证明。但是对于\(n\nmid k\)的情况就有点不知所措了。

不过别担心，这种情况也好做。注意到这是一个幂次逐渐递增的和式，恰好的是它的公比为\(\omega^{k}_{n}\)，于是运用求和公式我们有：

\[ \sum^{n-1}_{i=0}(\omega^{k}_{n})^i=\frac{1-(\omega^{k}_{n})^n}{1-\omega^{k}_{n}}=\frac{1-1}{1-\omega^{k}_{n}}=0 \]

这个引理在后面的快速插值运算中有着重要的作用。

有了这些东西，我们就可以开始进行优化了！

离散傅里叶/快速傅里叶变换（DFT/FFT）

我们要求的点值就是用这些单位复数根作为多项式的\(x\)代入，然后求出对应的\(f(x)\)值，这一操作称为离散傅里叶变换，快速傅里叶变换是对其进行加速。下面开始讲解FFT算法。

我们把多项式\(f(x)\)的奇偶项分开，构造如下的两个多项式的系数表示：
\[f_0=(a_0,a_2,\cdots,a_{n-2})\]
\[f_1=(a_1,a_3,\cdots,a_{n-1})\]
那么就有：
\[f(x)=f_0(x^2)+xf_1(x^2)\]
为什么要这么做呢？这就是根据单位复数根的性质来的。我们对于\(k<\frac{n}{2}\)，令\(x=\omega^k_n\)，那么有：
\[ f(\omega^k_n)=f_0(\omega^{2k}_n)+\omega^k_nf_1(\omega^{2k}_n)=f_0(\omega^{k}_\frac{n}{2})+\omega^k_nf_1(\omega^{k}_\frac{n}{2}) \]
\[ f(\omega^{k+\frac{n}{2}}_n)=f_0(\omega^{2k+n}_n)+\omega^{k+\frac{n}{2}}_nf_1(\omega^{2k+n}_n)=f_0(\omega^{k}_\frac{n}{2})-\omega^k_nf_1(\omega^{k}_\frac{n}{2}) \]
于是我们就可以递归求解了（不过我们一般用非递归做法求解，后面我会提到）。

这种求点值的方法的时间复杂度为\(O(n\log n)\)，从它的二分操作我们就容易得知了。

好了，我们到现在已经得到了点值，那么接下来怎么做呢？难道用拉格朗日插值法吗？当然不是，拉格朗日插值法现在用的话只会添麻烦，于是我们有一个更优秀的方法。

逆离散傅里叶变换（IDFT）

我们首先要构造几个矩阵，和这个多项式有一些关系，样子是这样的：

\[ \begin{pmatrix} 1 &\omega^0_n &\omega^0_n &\cdots &\omega^0_n \\ 1 &\omega^1_n &\omega^2_n &\cdots &\omega^{n-1}_n \\ 1 &\omega^2_n &\omega^4_n &\cdots &\omega^{2(n-1)}_n \\ \vdots &\vdots &\vdots &\ddots &\vdots &\\ 1 &\omega^{n-1}_n &\omega^{2(n-1)}_n &\cdots &\omega^{(n-1)^2}_n \\ \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n-1} \end{pmatrix} = \begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{n-1} \end{pmatrix} \]
注意到左边那个矩阵每行代表了一个\(x\)的取值，而每列又是逐一递增的幂次（其实就是每行都是一个相同的无系数多项式，只是\(x\)值不同罢了）。我们发现它与中间那个系数构成的矩阵相乘可以得到\(y\)值矩阵。

另外，像左边的那个矩阵的形式我们一般称为系数矩阵这里我们将其记作\(V_n\)。

为什么要构建这样的矩阵呢？因为我们要根据点值快速求系数。但是乍一看好像也是\(O(n^2)\)的算法，实际上并非如此，我们可以用上矩阵的性质将这个复杂的乘式变为我们熟悉的东西。

矩阵的性质：如果已知矩阵\(A\)和\(C\)且满足\(AB=C\),则\(B=A^{-1}C\)，其中\(A^{-1}\)表示矩阵\(A\)的逆。矩阵的逆与原矩阵相乘会得到单位矩阵。

按照这一性质，我们就可以很快地求出系数了，具体操作还在下面。

首先我们需要求出系数矩阵的逆。根据矩阵的性质我们可以得出：
\[ V_nV^{-1}_n=I_n= \begin{pmatrix} 1&&&&\\ &1&&&\\ &&1&&\\ &&&\ddots &\\ &&&&1 \end{pmatrix} \]
为了方便，我这里没写出\(0\)了，单位矩阵就只有主对角线有\(1\)，其它地方都是\(0\)。

可能你看到这个庞然大物的时候有点惊慌失措，不好怎么求矩阵的逆，但是不要担心，并不难，一步一步来你就明白了。

我们先只考虑主对角线的元素都变成\(1\)的情况，先对各个元素定好位再说，那么我们只需要使得逆矩阵的第\(i\)列与原矩阵的第\(i\)行的对应元素相乘后在求和等于\(1\)即可。那其它地方怎么办？求和引理为我们很好地解决了这个问题。

我们现在只需定位到逆矩阵的列，先不说得到\(1\)，按照求和引理我们先得到\(n\)。将逆矩阵第\(i\)列第\(j\)个元素与原矩阵第\(i\)行的第\(j\)个元素相乘后，其幂次\(k\)必须满足\(n\mid k\)，我们只需将\(k\)归零即可，最后乘得的那些\(n\)再除一下就得到了\(1\)。那其它地方的\(0\)呢，很简单，其它地方对应的行列是交错的，必然不会时时刻刻满足\(n\mid k\)，但是会满足幂次线性（线性+线性=线性）递增，导致其最终变成了\(0\)。

于是我们得到了：
\[ \frac{1}{n}\begin{pmatrix} 1 &\omega^0_n &\omega^0_n &\cdots &\omega^0_n \\ 1 &\omega^{-1}_n &\omega^{-2}_n &\cdots &\omega^{-(n-1)}_n \\ 1 &\omega^{-2}_n &\omega^{-4}_n &\cdots &\omega^{-2(n-1)}_n \\ \vdots &\vdots &\vdots &\ddots &\vdots &\\ 1 &\omega^{-(n-1)}_n &\omega^{-2(n-1)}_n &\cdots &\omega^{-(n-1)^2}_n \\ \end{pmatrix} \begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{n-1} \end{pmatrix} = \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n-1} \end{pmatrix} \]

都到这一步了，难道还要我们用平方暴力吗？不！

既然多项式可以变为系数矩阵，那么系数矩阵也可以变为多项式啊！我们只需要将FFT稍微修改一下，最后除以一个\(n\)就行了。真是妙不可言啊！

当然我们通常会使用非递归的方法来写FFT，我们称之为FFT的迭代实现，详见下方。

快速傅里叶变换的迭代实现

搬上经典的《算法导论》的图：

设\(n=2^k\)，我们注意到它们在第\(i\)层（根为第\(0\)层）的时候被按照其值的二进制第\(i\)位的异同被划分开来，划分开之后，它们位置的顺数第\(i\)个二进制位就被确定了，即第\(k-i-1\)位被确定。又因为划分到左侧的数第\(i\)位是\(0\)，划分之后位置的第\(k-i-1\)位也是\(0\)，右侧同理。于是它们最后被分配到的位置的二进制值与它们原来的二进制值是反过来的。

于是迭代法可以这么实现：我们就可以先求出它们反过来的二进制值，然后交换值，接着就按照FFT的算法进行求解。

而且，求反过来的二进制值也可以用递推的方法求解。

递推式就是：\(rev[i]=rev[i/2]/2\vee ((1\And i)/(l-2)),i< rev[i]\)

我们可以这么想：首先，\(i\)之前的\(rev\)都算过了，而\(rev[i/2]\)是\(i/2\)其反过来的二进制，那么除以\(2\)，就相当于原二进制乘以\(2\)。但是注意到乘以\(2\)后第\(0\)位必然是\(0\)，而\(i\)的第\(0\)位有可能是\(1\)，不过我们已经得到了这个数的反过来的二进制了（除了那一位不确定的），那我们只需要看一下这个\(i\)的第\(0\)位即可，然后将这个没倒过的位倒过去，或运算一下就得到了我们要的这个反过来的二进制，位置交换就好了。

交换好之后，我们就可以进行合并操作了，算好单位根，次数逐次递增即可。

至此，FFT教程也讲解完毕了。

代码实现：

#include <cstdio>
#include <cmath>
#include <algorithm>

using namespace std;

struct cn//虚数相关
{
    double a, b;
    cn(double a = 0, double b = 0) : a(a), b(b) {}
};
const double pi = 3.141592653589793;

cn operator + (const cn& a, const cn& b)
{
    return cn(a.a + b.a, a.b + b.b);
}

cn operator - (const cn& a, const cn& b)
{
    return cn(a.a - b.a, a.b - b.b);
}

cn operator * (const cn& a,const cn& b)
{
    return cn(a.a * b.a - a.b * b.b, a.a * b.b + a.b * b.a);
}

inline void dft(cn a[], int n, int l, int f)
{
    int rev[n];
    cn x, y;
    rev[0] = 0;
    for(register int i = 1; i < n; i += 1)
    {
        rev[i] = rev[i>>1] >> 1 | ((i&1) << l-1);//这里可以预处理，写在外面
        if(i < rev[i])
            swap(a[i], a[rev[i]]);
    }
    for(register int i = 1; i < n; i <<= 1)
    {
        cn wi(cos(pi / i), f * sin(pi / i));
        //注意这里本应是2pi/2i，因为2被约掉了
        //f=-1的时候只要次数变负即可，再根据欧拉公式得到只要sin前面加个负号即可
        for(register int j = 0; j < n; j += i<<1)//正因为是2n，所以要乘以2
        {
            cn w(1, 0);
            for(register int k = 0; k < i; k += 1)//这里不需要变i，因为k本来就小于deg的一半
            {
                x = a[j+k];
                y = a[j+k+i];
                a[j+k] = x + w * y;
                a[j+k+i] = x - w * y;
                w = w * wi;//旋转
            }
        }
    }
    if(f == -1)
        for(register int i = 0; i < n; i += 1)
            a[i].a /= n, a[i].b /= n;
}

int main()
{
    int n, m, N = 1, l = 0;
    double x;
    scanf("%d%d", &n, &m);
    int len = n+m+1;
    for(; N < len; N <<= 1)
        l += 1;
    cn a[N], b[N];
    for(register int i = 0; i <= n; i += 1)
    {
        scanf("%lf", &x);
        a[i] = cn(x, 0);
    }
    for(register int i = 0; i <= m; i += 1)
    {
        scanf("%lf", &x);
        b[i] = cn(x, 0);
    }
    dft(a, N, l, 1);
    dft(b, N, l, 1);
    //系数进去，y值出来
    //因为递归到底层就没有x了
    cn c[N];
    for(register int i = 0; i < N; i += 1)
        c[i] = a[i] * b[i];
    dft(c, N, l, -1);
    //求逆可以利用差不多的方法
    for(register int i = 0; i < len; i += 1)
        printf("%d ", (int)(c[i].a + 0.5));
    //注意精度
    return 0;
}

习题推荐

2017湖南/安徽省选题:礼物

假设旋转后的数列为：\({a},{b}\)，第一个手环亮度调整为\(c\)，那么按题目要求就是使得如下式子的值最小：

\[\sum^{n}_{i=0}(a_i -b_i +c)^2\]

之后我们再将里面的式子拆开，得（这应该初中就学过了公式了吧）：

\[\sum^{n}_{i=0}(a^2_i+b^2_i+c^2-2a_ib_i+2a_ic-2b_ic)\]

然后将求和分开，整理得：

\[\sum^{n}_{i=0}a_i^2 +\sum^{n}_{i=0}b_i^2 -2\sum^{n}_{i=0}(a_ib_i)+nc^2+2(\sum^{n}_{i=0}a_i-\sum^{n}_{i=0}b_i)c\]

可以发现，关于调整亮度的\(c\)仅仅是一个二次函数，其系数是确定的，要使和式最小，只需求出二次函数的最小值。并且注意到最前面两个平方和也是确定的，那么不能直接确定的就是。中间的那个乘积之和。

我们再把这个乘积之和提出来看，我们只需保证它最大就好了，我们将\(b\)环扩展到\(2n-1\)的长度（也就是把后面循环写出来），以便我们计算不同的旋转方案。假设我们第二个环从第\(k\)个开始，那么我们得到：

\[\sum^{n}_{i=0}(a_ib_{i+k-1})\]

但是这并不能用多项式加速，因为它不满足卷积的形式——两乘积的下标和为定值。

于是有一种方法就是将\(a\)的系数反过来，这种操作真令人佩服，即：

\[\sum^{n}_{i=0}(a_{n-i}b_{i+k-1})\]

虽然反转了系数，但是乘的是同一个东西，这是没问题的，因此我们就可以在复杂度为\(O(n\log n)\)的时间内完成计算，然后在线性查找最大值即可。

算完这些分步之后，将答案加加减减就可以了。

怎么样，多项式有意思吧！

然后，还有一些习题老师推荐的我也推荐给大家，尽管我还没做，但确实是好题。

Luogu P3338:力，好像是电场力的公式诶

Luogu P4199:万径人踪灭，字符串+多项式综合题

后置知识

本文提到的一些高深数学知识，列举如下：

弧度制（单位：\(rad\)）：\(\pi =180\)°

定义：圆上\(1rad\)角所对的弧长等于该圆的半径的长。

又因为\(C=2\pi r\)，故\(\pi =180\)°。

虚数单位：\(i=\sqrt{-1}\)

理解：我们可以考虑一个实数轴，对于每一个实数\(x\)与\(-1\)相乘的运算，我们都看作是\(x\)绕原点旋转了\(\pi\)，那么旋转\(2\pi\)就是\(-1\times(-1)=1\)。那么旋转\(\frac{\pi}{2}\)呢？于是就产生了虚数单位，它与实数组成复数集。它是跨出实数轴的一个东西，不存在于实数中，我们可以把复数集看作是一个二维平面的东西（我们称之为复平面），注意，仅“看作”。

自然对数的底数（自然常数）：\(e\)

其中有一个定义的求值公式为：
\[e=\sum^{\infty}_{i=0}\frac{1}{i!}\approx 2.71828\]
注意\(0!=1\)。这个值没必要记，了解就好了，有兴趣的同学可以了解一下它的指母函数和正余弦的指母函数，对理解欧拉公式有帮助（虽然指母函数也有些抽象）。

矩阵的话，在文章里面讲得很清楚了，这里就不写了吧2333（其实是博主懒）。

数学中的与或非：

与：\(\wedge\)

或：\(\vee\)

非：\(\urcorner\)