FFT • CoderXL's Blog

使用复数（单位根）和分治（蝶形运算）加速多项式乘法。

问题#

计算两个次数为 $n$ 和 $m$ 的多项式 $A(x),B(x)$ 的乘积，结果为一个 $n+m$ 次多项式 $C(x)$ .

朴素#

朴素算法使用系数对系数乘法，要算出结果多项式中的每一项的系数都需要 $O(n)$ ，整体复杂度 $O(n^2)$ .

优化#

考虑优化。众所周知，分治能够优化很多问题，而其适用条件是：父问题可以拆成相似的几个子问题，且当子问题的答案已经求出时，能够快速地求解（或称“合并”）出父问题的答案。

乍一看，多项式乘法中似乎没有什么方法，能够施展“先拆分后合并”的思想。但是下面这种方式证明了分治是可行的：

Karatsuba 分治#

可以达到 $O(n^{\log_2 3})\approx O(n^{1.585})$ 的复杂度。

流程：

先拆分：将 $n$ 次多项式 $P(x)$ 按照次数从低到高，以 $x^w$ 为界拆成两半： $P(x)=P_0(x)+x^w P_1(x)$ ，而两个子多项式分别是 $w-1$ 和 $n-w$ 次的. 这个问题中，我们将待乘的两个多项式均如是拆分，得到 $A_0(x),A_1(x),B_0(x),B_1(x)$ .
再合并：朴素的算法流程是 $C=A_0B_0+x^w A_1B_0+x^w A_0B_1+x^{2w} A_1B_1$ ，需要计算 $4$ 次子多项式乘法。使用 Master Theorem 分析得知其复杂度仍为 $O(n^2)$ . 但是，我们可以利用一个数学技巧： $A_1B_0 + A_0B_1=(A_0+A_1)(B_0+B_1)-A_0B_0 - A_1B_1$ ，把中间的两次多项式乘法变成两次多项式加法、一次多项式乘法、两次多项式减法，复杂度来到 $O(n^{\log_2 3})$ .

注意到 $w$ 的选取可能影响到算法的优劣，当 $A(x),B(x)$ 次数相同时，取 $\lceil {n\over2} \rceil$ 是个不错的选择；但当二者长度不相同时，应当取 $\lceil {\max\{n,m\} \over2} \rceil$ ，但此时对于次数较小的一侧会有计算上的浪费。

更多、更快的类似分治优化也存在，但是实现复杂度较高，且其优化效果不如即将介绍的这个全新思路。

FFT#

上述优化是对系数域进行的，离不开系数之间的卷积。但如果利用拉格朗日插值法的思想，唯一确定一个 $n+m$ 多项式只需要 $n+m+1$ 个点即可，而每个点只需要在原来的两个多项式中采样后相乘，复杂度为 $O(n)$ .

这启发我们将系数乘法转化为点值乘法，完成运算后再逆转换回去，或许能够大幅加快算法速度。

然而，要获得每个点的点值，我们会考虑直接代入 $n+m+1$ 个点，但是每个点的函数值计算复杂度仍然是 $O(n)$ ，因此由多项式得到点值这一步本身就是 $O(n^2)$ 的。但是，这一步通过一些技巧能够大幅优化。

单位根#

在复平面上，单位圆上的点有很多特殊的性质。

而所谓 $n$ 次单位根 $z_n$ 是指满足 $z_n^n=1$ 的点。不难发现有 $n$ 个这样的点且他们均在单位圆上，其表达式是 $z_{n,k}=e^{{2\pi k\over n}i},\ k\in\{0,1,\dots,n-1\}$ . 而由于所有的 $z_{n,k}$ 都能用 $z_{n,1}^k$ 表示，我们又将 $z_{n,1}$ 简记为 $z_n$ ，进而所有 $n$ 次单位根都可以记为 $z_n^k,\ k\in\{0,1,\dots,n-1\}$ .

回到上文，我们希望取（至少） $n+m+1$ 个不同的点并计算多项式在这些点处的值。如果我们令 $p=2^{\lceil \log_2(m+n+1) \rceil}$ ，然后取 $p$ 次单位根，那么总共可以得到 $p\ge n+m+1$ 个点（这里 $p=\Theta(n)$ ，即 $p$ 和 $n$ 是同规模的）。而我们将看到，计算多项式在这些点的值的时候，由于单位根特殊的数学性质，很多数值可以复用，因此能够大幅简化计算。

具体地，计算一个具有 $2n$ 项的多项式 $f(x)=a_0+a_1x+\cdots+a_{2n-1}x^{2n-1}$ 在 $x=z_p^k$ 处的值时，我们可以将原式按次数的奇偶性拆分：

\begin{aligned} f(x)&=(a_0+a_2x^2+\cdots+a_{2n-2}x^{2n-2})+(a_1x+a_3x^3+\cdots+a_{2n-1}x^{2n-1})\\ &=(a_0+a_2x^2+\cdots+a_{2n-2}x^{2n-2})+x(a_1+a_3x^2+\cdots+a_{2n-1}x^{2n-2}) \end{aligned}

然后令 $g(x)=a_0+a_2x+\cdots+a_{2n-2}x^{n-1}$ ， $h(x)=a_1+a_3x+\cdots+a_{2n-1}x^{n-1}$ ，得到 $f(x)=g(x^2)+xh(x^2)$ .

此时，我们发现 $g(x^2)$ 和 $h(x^2)$ 这两个需要计算的部分竟然是偶函数，也就意味着，一旦计算出了 $g(x^2)$ 和 $h(x^2)$ ，不但可以得到 $f(x)$ ，还可以把 $f(-x)=g(x^2)-xh(x^2)$ 一并算出。而 $-z_p^k=z_p^{k+{p/2}}$ 正是我们需要计算的另一个单位根。因此，我们利用了单位根的性质，成功做到了算一遍得两个取值。

注意到 $(z_p^k)^2=z_p^{2k}$ 仍然是单位根，因此 $g(x^2)$ 和 $h(x^2)$ 的计算可以很方便地递归进行下去，但这要求 $f(x)$ 的项数是 $2$ 的幂。一个常见的做法是给 $f(x)$ 的高次项补 $0$ 到最低的 $2^h-1$ 次。

上述做法相当于，将「计算 $2n$ 项多项式在 $p$ 个点处的函数值」这个任务，转化为了 $2$ 个「计算 $n$ 项多项式在 ${p\over 2}$ 个点处的函数值」的任务，然后进行复杂度为 $O(n)$ 的合并。总复杂度 $O(n\log n)$ .

刚刚所描述的，就是快速离散傅里叶变换（FFT）。

逐点相乘#

通过上述方法， $O(n\log n)$ 地得到了 $p$ 个点处两个多项式的函数值之后，可以进行逐点对应相乘，复杂度 $O(n)$ .

从点值还原为系数#

我们刚刚利用 FFT 将多项式从系数形式变换为点值形式，随后逐点相乘得到了

C(z_p^k)=A(z_p^k)\cdot B(z_p^k),\quad k\in{0,1,\dots,p-1}

而我们真正想要的，是 $C(x)$ 的系数形式。

由于 $p\ge n+m+1$ ，并且 $C(x)$ 的次数小于等于 $n+m$ ，因此 $C(x)$ 可以被这 $p$ 个点唯一确定。换言之，我们只需做一次插值即可还原出 $C(x)$ 的系数。

但朴素的拉格朗日插值仍然是 $O(n^2)$ 的，因此我们必须进一步利用单位根的性质。

事实上，FFT 所做的事情等价于离散傅里叶变换（DFT）：

Y_k=\sum_{j=0}^{p-1} a_j z_p^{kj}

而逆变换（IDFT）为：

a_j={1\over p}\sum_{k=0}^{p-1} Y_k z_p^{-kj}

注意到二者的形式几乎完全一致，只是指数取了负号，并多了一个系数 ${1\over p}$ .

因此，我们可以用几乎同样的分治结构完成逆变换：只需在递归合并时，将旋转因子 $z_p^k$ 换成其逆元 $z_p^{-k}$ ，最后将所有系数统一除以 $p$ 即可。

这个过程被称为逆 FFT（IFFT），其复杂度同样为 $O(n\log n)$ .

至此，多项式乘法的完整流程为：

选取 $p=2^{\lceil \log_2(m+n+1) \rceil}$ 并对两多项式补零至 $p$ 项；
对 $A,B$ 分别进行 FFT 得到点值表示；
逐点相乘得到 $C$ 的点值表示；
对 $C$ 进行 IFFT 得到系数表示（并对结果除以 $p$ ）；

整体复杂度为 $O(n\log n)$ ，相比朴素的 $O(n^2)$ 有数量级的提升。

最后需要注意，FFT 使用浮点数复数运算，会引入误差，因此输出系数通常需要进行四舍五入处理。对于系数较大的情况，有时还需将系数拆分为高低位分别卷积以避免精度爆炸。

蝶形优化#

可以省去递归分治时使用的 swap 数组，减小空间复杂度，具体可以参考模板代码。

模板代码#

//
//  main.cpp
//  P3803
//
//  Created by Leo Xia on 2026/3/29.
//

#include <bits/stdc++.h>
using namespace std;

const int N = 1000010<<3;//空间复杂度可能达到 (n+m)*4 ≈ 8n
const double Pi = acos(-1.0);
struct C//Complex number
{
    double x,y;
    C()=default;
    C(double x_,double y_):x(x_),y(y_){}
    C operator+(const C&b)const
    {return {x+b.x,y+b.y};}
    C operator-(const C&b)const
    {return {x-b.x,y-b.y};}
    C operator*(const C&b)const
    {return {x*b.x-y*b.y,x*b.y+y*b.x};}
}a[N],b[N];

int n,m;
int lim=1,l;
int r[N];

void fft(C p[], int sgn)
{
    for(int i=0;i<lim;i++)
        if(i<r[i])swap(p[i],p[r[i]]);
    for(int mid=1;mid<lim;mid<<=1)
    {
        C wn(cos(Pi/mid),sgn*sin(Pi/mid));
        for(int len=mid<<1,j=0;j<lim;j+=len)
        {
            C w(1,0);
            for(int k=0;k<mid;k++,w=w*wn)
            {
                C x=p[j+k],y=w*p[j+k+mid];
                p[j+k]=x+y;
                p[j+k+mid]=x-y;
            }
        }
    }
}

int main()
{
    cin>>n>>m;
    for(int i=0;i<=n;i++)cin>>a[i].x;
    for(int i=0;i<=m;i++)cin>>b[i].x;
    while(lim<=n+m)
    {
        lim<<=1;
        l++;
    }
    for(int i=0;i<lim;i++)
    {
        r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    }
    fft(a,1);
    fft(b,1);
    for(int i=0;i<lim;i++)
        a[i]=a[i]*b[i];
    fft(a,-1);
    for(int i=0;i<=n+m;i++)
        cout<<(int)(a[i].x / lim + 0.5)<<' ';//误差处理
    return 0;
}

cpp