多项式卷积

给定两个至多 $n-1$ 次的多项式 $$A(x)=\sum _{k=0}^{n-1}{a}_k{x}^k,B(x)=\sum _{k=0}^{n-1}{b}_k{x}^k$$ 如何快速计算两者的卷积，即它们相乘得到的多项式的系数？ $$\begin{array}{rl}(A\ast B)(x)=A(x)B(x)& =\sum _{i=0}^{n-1}{a}_i{x}^i\sum _{j=0}^{n-1}{b}_j{x}^j\\ & =\sum _{k=0}^{2n-2}{x}^k\sum _{i+j=k}{a}_i{b}_j\end{array}$$ 上式给出了 $O(n^2)$ 的朴素做法．

系数 - 点值 - 系数

• 代入任意 $x$ 可得到多项式在 $x$ 处的点值

• 点值意义下的多项式乘法是 $O(n)$ 的 $$(A\ast B)(x)=A(x)B(x)$$

• $n$ 点确定一个至多 $n-1$ 次的多项式¹

若计算至多 $n-1$ 次的多项式的某 $n$ 个点值存在快速算法，而通过多项式的某 $n$ 个点值确定原多项式系数亦存在快速算法，就有通过 $O(n)$ 的点值乘法加速多项式乘法计算的可能．

系数 - 点值 - 系数 - 快速转换？

• 朴素计算任意指定 $n$ 个位置点值需要 $O(n^2)$．

• Lagrange[1] 插值给出了 $O(n^2)$ 将任意位置 $n$ 个点值还原为多项式系数的算法．

• 能否选取 $n$ 个特殊的点值使系数 - 点值、点值 - 系数的变换支持快速计算？

Discrete Fourier Transform

离散傅里叶变换（Discrete Fourier Transform, DFT）接受一个至多 $n-1$ 次的多项式的 $n$ 个系数，将复数域上的 $n$ 个 $n$ 次单位根代入系数表达式以得到给定多项式的 $n$ 个点值． $$a_0,a_1,\dots ,a_{n-1}⟶A(1),A({\omega }_n),\dots ,A({\omega }_n^{n-1})$$

得益于单位根的特殊运算性质，二者均有被称为快速傅里叶变换（Fast Fourier Transform, FFT）的快速算法．

复数域单位根

复数域上的 $n$ 个 $n$ 次单位根² $${\omega }_n^k:=e^{\frac{2\pi k}{n}i}=\cos \frac{2\pi k}{n}+i\sin \frac{2\pi k}{n},k\in [n]$$ 是复平面单位圆上的 $n$ 等分点，易验证它们是复数域中唯一满足方程 $z^n=1$ 的解．

所有单位根模长均为 $1$．第 $k$ 个单位根的辐角为 $\frac{2\pi k}{n}$．复数乘法”模长相乘，辐角相加”的性质告诉我们 ${\omega }_n^i{\omega }_n^j={\omega }_n^{i+j}$．

复数域单位根 - 三个重要性质³

消去 / 折半引理将在 FFT 的推导中使用．

复数域单位根 - 三个重要性质

求和引理的证明使用了等比数列求和公式．将在 IDFT 的推导中用到．

Fast Fourier Transform

考虑将至多 $2n-1$ 次的待变换多项式 $A(x)=\sum _{k=0}^{2n-1}{a}_k{x}^k$ 奇偶分项两个至多 $n-1$ 次的多项式 $A(x)=A_0(x^2)+x{A}_1(x^2)$，其中 $$A_0(x)=\sum _{k=0}^{n-1}{a}_{2k}{x}^k,A_1(x)=\sum _{k=0}^{n-1}{a}_{2k+1}{x}^k$$ 代入 $x={\omega }_{2n}^k,k\in [2n]$，用单位根消去或折半引理（${\omega }_{2n}^{2k}={\omega }_n^k$）得 $$A({\omega }_{2n}^k)=A_0({\omega }_n^k)+{\omega }_{2n}^k{A}_1({\omega }_n^k)$$ 再用 ${\omega }_{2n}^n=-1$ 得 $$\begin{array}{rl}A({\omega }_{2n}^k)& =A_0({\omega }_n^k)+{\omega }_{2n}^k{A}_1({\omega }_n^k)\\ A({\omega }_{2n}^{n+k})& =A_0({\omega }_n^k)-{\omega }_{2n}^k{A}_1({\omega }_n^k)\end{array},k\in [n]$$ 原多项式 $A(x)$ 规模为 $2n$ 的 DFT 转化为规模为 $n$ 的 $A_0(x)$ 和 $A_1(x)$ 的 DFT．递归计算就可得到 $O(n\log n)$ 的算法．

DFT 的矩阵表示

记 $$\begin{array}{c}\mathit{a}=\left(\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{n-1}\end{array}\right),F={\left({\omega }_n^{ij}\right)}_{(i,j)\in n\times n}=\left(\begin{array}{cccc}1& 1& \dots & 1\\ 1& {\omega }_n& \dots & {\omega }_n^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{n-1}& \dots & {\omega }_n^{(n-1)(n-1)}\end{array}\right)\end{array}$$ 则 DFT 的变换结果（给定系数 $a_0,a_1,\dots ,a_{n-1}$ 的多项式在 $n$ 个 $n$ 次单位根处的点值）可表示为 $$\hat{\mathit{a}}=\left(\begin{array}{c}A(1)\\ A({\omega }_n)\\ ⋮\\ A({\omega }_n^{n-1})\end{array}\right)=F\mathit{a}$$

IDFT

由单位根的消去引理可证，DFT 矩阵 $F$ 的逆矩阵为⁴ $$F^{-1}=\frac{1}{n}{\left({\omega }_n^{-ij}\right)}_{(i,j)\in n\times n}=\frac{1}{n}\left(\begin{array}{cccc}1& 1& \dots & 1\\ 1& {\omega }_n^{-1}& \dots & {\omega }_n^{-(n-1)}\\ ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{-(n-1)}& \dots & {\omega }_n^{-(n-1)(n-1)}\end{array}\right)$$ 于是 IDFT 的过程可表示为 $$\mathit{a}=F^{-1}\hat{\mathit{a}}$$

故快速计算 IDFT 的方法与 FFT 几乎一致，只需将计算 DFT 时使用的本原单位根 ${\omega }_n$ 替换为 ${\omega }_n^{-1}$ 并对最终结果除以 $n$ 即可．

梳理

• 怎么计算卷积？

  • 把至多 $n-1$ 次的多项式 $A(x)$ 和至多 $m-1$ 次的多项式 $B(x)$ 写成至多 $n+m-2$ 次的多项式（高位补 $0$）．为计算 FFT 方便，还要继续补 $0$ 至一个大于其次数的 $2$ 的幂．
  • 对 $A(x)$ 和 $B(x)$ 多点求值．
  • 把两个多项式的点值逐点相乘．
  • 多点插值还原 $(A\ast B)(x)$ 的系数．

• 怎么快速求值？

  • 选点选单位根就是 DFT．
  • 消去引理和折半引理使我们可以用 FFT 算法递归地计算 DFT．
  • 推导已经给出了递归的写法，之后还会介绍常数更优的迭代实现．

• 怎么快速插值？

  • 求和引理给出了 DFT 矩阵的逆矩阵．
  • 计算方法很相似，最后逐项除掉一个规模．

FFT 递归实现 - DFT 部分

方便起见，我们只处理 $n$ 为 $2$ 的幂的情形．以下 C 风格的代码实现了递归的 DFT 和 IDFT．

#include<bits/stdc++.h>
#include<complex>
using namespace std;
typedef long long ll; typedef complex<double> CP;
const ll MXN=4E6+5; const double PI=3.14159265358979323846l;
CP tmp[MXN];
void _DFT(CP A[],ll n,ll typ){
    n/=2; if(n==0) return;
    for(ll k=0;k<n;k++) tmp[k]=A[2*k],tmp[n+k]=A[2*k+1];
    for(ll k=0;k<2*n;k++) A[k]=tmp[k];
    _DFT(A,n,typ); _DFT(A+n,n,typ);
    CP w(cos(2*PI/(2*n)),typ*sin(2*PI/(2*n))), wk=1;
    for(ll k=0;k<n;k++){
        tmp[  k]=A[k]+wk*A[n+k];
        tmp[n+k]=A[k]-wk*A[n+k];
        wk*=w;
    } for(ll k=0;k<2*n;k++) A[k]=tmp[k];
} void DFT(CP A[],ll n,ll typ){
    _DFT(A,n,typ); if(typ==-1) for(ll i=0;i<n;i++) A[i]*=1.0l/n;
}

FFT 递归实现 - 卷积部分

// alternatively, use std::__lg() in GCC
ll log2ceil(ll n){ll cnt=0; for(ll t=1;t<n;t<<=1) cnt++; return cnt;} 
CP A[MXN],B[MXN],C[MXN]; ll outC[MXN];
ll* conv(ll inA[],ll aN,ll inB[],ll bN){
    ll n=1LL<<log2ceil(aN+bN-1);
    for(ll i=0;i<aN;i++) A[i]=inA[i];
    for(ll i=0;i<bN;i++) B[i]=inB[i];
    DFT(A,n,1); DFT(B,n,1);
    for(ll i=0;i<n;i++) C[i]=A[i]*B[i];
    DFT(C,n,-1); for(ll i=0;i<n;i++) outC[i]=round(C[i].real());
    return outC;
}

Drawbacks?

• 递归实现慢
• 临时数组丑
• 封装性为零

FFT 迭代

迭代地实现 FFT 不仅在常数上更加优秀，亦更便于使用 C++ 的容器进行封装．这并不困难，只需自底向上模拟 FFT 递归过程即可．

唯一的问题——最底层的顺序？

来观察一轮 $2^3$-FFT 自顶向下的置换过程 $$\begin{array}{ccccccccc}{2}^3& 0/000& 1/001& 2/010& 3/011& 4/100& 5/101& 6/110& 7/111\\ 2^2& 0/000& 2/010& 4/100& 6/110& 1/001& 3/011& 5/101& 7/111\\ 2^1& 0/000& 4/100& 2/010& 6/110& 1/001& 5/101& 3/011& 7/111\\ 2^0& 0/000& 4/100& 2/010& 6/110& 1/001& 5/101& 3/011& 7/111\end{array}$$

你发现了什么？

蝶形运算

在计算点值前，$2^n$-FFT 事实上完成了一次 $n$-位逆序置换．分解来看，规模为 $2^k$ 的层的置换完成了对 $2^{n-k}$ 对应二进制位的分类．

我们有 $O(n)$ 的递推方法获得这一置换．

void spawnrev(ll n){
    rev[0]=0;
    for(ll i=1;i<(1<<n);i++)
        rev[i]=(rev[i>>1]>>1)+((i&1)<<(n-1));
}

FFT 迭代实现

void DFT(CP A[],ll n,ll typ){ // rev[] should be spawned in advance
    for(ll i=0;i<n;i++) if(i<rev[i]) swap(A[i],A[rev[i]]); // a one-to-one permutation
    for(ll hf=1;hf<n;hf*=2){
        CP w(cos(2*PI/(2*n)),typ*sin(2*PI/(2*n))), wk=1;
        for(ll i=0;i<n;i+=hf*2){
            CP wk=1;
            for(ll k=0;k<hf;k++){
                CP x=A[i+k],y=wk*A[i+hf+k];
                A[i+k]=x+y; A[i+hf+k]=x-y;
                wk=wk*w;
            }
        }
    }
    if(typ==-1) for(ll i=0;i<n;i++) A[i]*=1.0l/n;
}

请自行实现更易用的容器封装版本．

NTT 速成⁵

FFT 的缺点？浮点数带来的大常数与精度问题．

$998244353=2^{23}\times 7\times 17+1$，同时是一个质数．

我们指出，在系数和点值模 $p=998244353$ 的意义下，当规模 $n\le p$ 时⁶，至多 $n-1$ 次的多项式仍可由其 $n$ 个点值唯一确定，故仍可使用系数-点值-系数的方法求得多项式卷积．

我们指出，对于满足 $n\mid p-1$ 的 $n$，依 ${\omega }_n:=3^{\frac{p-1}{n}}$ 定义的 ${\omega }_n$ 在模 $p$ 意义下与复数域中定义的 ${\omega }_n$ 发挥相同的作用，仍可进行规模至多为 $2^{23}$ 的 FFT 作为模 $p=998244353$ 意义下 NTT 的快速算法．

只需修改单位根定义，把复数运算改为整数取模，就得到了能算 NTT 的 FFT 的实现．

const PR=3,MOD=998244353;
ll w=qpow(PR,(MOD-1)/(hf*2)); if(typ==-1) w=inv(w);

FFT/NTT in a nutshell - 小结：概念区分

• 关于 DFT

  • Discrete Fourier Transform, DFT, 离散傅里叶变换
  • Fast Fourier Transform, FFT, 快速傅里叶变换
  • FFT 是计算 DFT 的快速算法

• 关于 NTT

  • Number Theoretic Transform, NTT, 数论变换
  • FFT 在复数域上的多项式环 $\mathbb{C}[x]$ 中进行，而 NTT 在模 $p$ 剩余类域 ${\mathbb{Z}}_p[x]$ 上进行
  • 快速计算 DFT / IDFT, NTT 都用 FFT，故一般不使用”FNTT”的说法

基本应用

大整数乘法

十进制数可拆解为多项式表示，计算卷积后处理进位即可．由于数字最大只是 $9$，合理数据范围下卷起来不会爆模数，直接用 NTT 实现即可．

背包计数

两个背包的合并就是多项式卷积．

滤波器

反转多项式的系数数组再做卷积，可以快速得到两个多项式滑动窗口式的内积．

位运算

有些位运算可以写成卷积的形式．模 $2$ 意义下”异或”是加法，“与”是乘法，“或”可以通过取反转化为”与”．

字符串

通过巧妙设定字符串距离函数，FFT 可解决更广泛的字符串匹配问题．

值域次数化

当值域较小时，将待计算的值放在多项式次数上统计贡献次数，可以绕开某些极难求解的数值问题．例如 Vandermonde 行列式的快速计算．

基本应用 - 分治 FFT

对多个长度相同的多项式的卷积，分治地卷起来可降低时间复杂度．长度不一时，挑小的先卷也可减小常数（用堆维护）．

另有一种 CDQ 风格的分治 FFT．CDQ 长于处理带偏序的二元点对贡献，在处理形如 $c_k=\sum _{i>j}[i+j=k]a_i{b}_j$ 的带偏序卷积时可以应用．

当卷积的前后项存在依赖关系时，也可使用此法保证处理顺序恰当．然而此类依赖问题往往也可通过解生成函数方程的方法求得封闭形式．

时间复杂度均为 $$T(n)=2T\left(\frac{n}{2}\right)+O(n\log n)=O(n{\log }^{2}n)$$

生成函数初步 - 导言

生成函数是一种对数列的操作技巧．通过将数列表示为多项式或形式幂级数，数列间复杂的和式操作可用简单函数的乘法、复合等运算进行表示，从而大大降低了数列变换技巧的使用门槛．

生成函数在组合数学中应用广泛，且生成函数的部分操作在组合意义下也有较为直观的理解．本节将带大家初窥其中的奥妙．限于篇幅和主讲人能力，我们仅以题带点地讲解，期冀为大家建立构造生成函数的直觉．请感兴趣的同学下来做进一步研究．

熟悉 Taylor 展开的同学或能较快上手此部分内容．

Ordinary Generating Function

序列 $a_n$ 的普通生成函数（Ordinary Generating Function, OGF）定义为其对应的多项式（形式幂级数）$A(x)=\sum _{n=0}^{+\mathrm{\infty }}{a}_n{x}^n$．这样书写的目的是为了便于进行数列间的卷积操作．特别的，Taylor 展开（或广义二项式定理）风格的 $\frac{1}{1-x}=1+x+x^2+\dots$ 也是常用技巧之一⁷．

OGF 组合意义

OGF 相乘，是背包，是卷积，是两块无标号组合对象的有序拼接．

Exponential Generating Function

序列 $a_n$ 的指数生成函数（Exponential Generating Function, EGF）定义为 $$A(x)=\sum _{n=0}^{+\mathrm{\infty }}{a}_n\frac{x^n}{n!}$$

EGF 组合意义

观察两个 EGF 的乘积 $$\begin{array}{rl}A(x)B(x)& =\sum _{i=0}^{+\mathrm{\infty }}{a}_i\frac{x^i}{i!}\sum _{j=0}^{+\mathrm{\infty }}{b}_j\frac{x^j}{j!}\\ & =\sum _{i=0}^{+\mathrm{\infty }}\sum _{j=0}^{+\mathrm{\infty }}{a}_i{b}_j\frac{x^{i+j}}{i!j!}\\ & =\sum _{n=0}^{+\mathrm{\infty }}\sum _{i=0}^n{a}_i{b}_{n-i}\frac{x^n}{i!(n-i)!}\\ & =\sum _{n=0}^{+\mathrm{\infty }}\sum _{i=0}^n{a}_i{b}_{n-i}(\genfrac{}{}{0ex}{}{n}{i})\frac{x^n}{n!}\end{array}$$ EGF 的乘积，相当于两块有标号组合对象的有序拼接．

EGF 组合意义

$$A(x)B(x)=\sum _{i=0}^{+\mathrm{\infty }}{a}_i\frac{x^i}{i!}\sum _{j=0}^{+\mathrm{\infty }}{b}_j\frac{x^j}{j!}=\sum _{n=0}^{+\mathrm{\infty }}\sum _{i=0}^n{a}_i{b}_{n-i}(\genfrac{}{}{0ex}{}{n}{i})\frac{x^n}{n!}$$

NTT 原理 - 导言

本节介绍 NTT 的原理．

FFT 加速卷积算法的核心，一是多项式的求值插值原理，二是单位根带来的分治快速算法．我们将在本节中证明，模 $p$ 剩余类域 ${\mathbb{Z}}_p$ 中，只要 $n\le p$，仍可通过 $n$ 个点值唯一确定至多 $n-1$ 次的 ${\mathbb{Z}}_p$ 上的多项式，这保证了系数-点值-系数方法的正确性；此外，基于原根构造出的 ${\mathbb{Z}}_p$ 上的本原单位根，使得使用与 FFT 相同的分治算法加速 NTT 的计算成为可能．

NTT 原理涉及原根等数论内容．本讲的目标是建立 DFT 变换和 FFT 算法的通用数学框架，而非具体研究其某一特例．故我们只讲解 NTT 所需的基础数论知识，无关的细节则略过处理．对数论感兴趣的同学可前往 OI Wiki 学习．

模 $p$ 剩余类域 ${\mathbb{Z}}_p$

模 $p$ 剩余类域 ${\mathbb{Z}}_p$ 是刻画取模运算下整数加法、乘法运算规律的代数结构，其良定义性由 $$\{\begin{array}{l}{a}_1\equiv a_2(\mathrm{mod}p)\\ b_1\equiv b_2(\mathrm{mod}p)\end{array}⟹\{\begin{array}{l}{a}_1+b_1\equiv a_2+b_2(\mathrm{mod}p)\\ a_1{b}_1\equiv a_2{b}_2(\mathrm{mod}p)\end{array}$$ 保证．

模 $p$ 剩余类域的最大特点是其中每个元素都存在（唯一）逆元．这一点是数论中 Bézout 定理的直接应用．

模 $p$ 剩余类域的另一优良性质是其满足消去律，即 $$ab\equiv 0(\mathrm{mod}p)⟹a\equiv 0\vee b\equiv 0(\mathrm{mod}p)$$ 而模一般的数 $m$ 形成的代数结构则不满足此性质．

${\mathbb{Z}}_p[x]$ 上的多点插值

之后记系数均在 ${\mathbb{Z}}_p$ 上的多项式组成的集合为 ${\mathbb{Z}}_p[x]$．

${\mathbb{Z}}_p[x]$ 上的多点插值

推论告诉我们，欲确定 ${\mathbb{Z}}_p$ 上的一个至多 $n-1$ 次的多项式，只需找到其 $n$ 个模 $p$ 意义下的点值即可．因此，使用系数 - 点值 - 系数方法加速 ${\mathbb{Z}}_p[x]$ 上的卷积是可行的．

本原单位根

回顾 DFT 中复数域 $\mathbb{C}$ 上单位根的定义 $${\omega }_n^k:=e^{\frac{2\pi k}{n}i}=\cos \frac{2\pi k}{n}+i\sin \frac{2\pi k}{n},k\in [n]$$ 其中 ${\omega }_n$ 是最重要的单位根，因为它生成了所有其它的单位根．

事实上所有的 ${\omega }_n^k:gcd(n,k)=1$ 也都是重要的，它们的 $0$ 至 $n-1$ 次方也能生成所有的单位根．这也是数论中 Bézout 定理的直接应用．

我们把这一类重要的单位根称为本原单位根．抽象的来说，$n$ 次本原单位根 ${\omega }_n$ 可定义如下 $${\omega }_n^n=1;{\omega }_n^k\ne 1,k=1,2,\dots ,n-1$$

${\mathbb{Z}}_m$ 上的阶

如何定义 ${\mathbb{Z}}_m$ 上的“本原单位根”呢？

称在模 $m$ 意义下使得 $a^n=1$ 的最小正整数 $n$ 为 $a$ 的阶，记为 ${\mathrm{ord}}_m(a)$．在不至混淆的情况下，也简记为 $\mathrm{ord}(a)$．

• 阶最高有多高？

Euler 定理的证明 - 简化剩余系

${\mathbb{Z}}_m^R:=\{a\in {\mathbb{Z}}_m:gcd(a,m)=1\}$ 被称为 ${\mathbb{Z}}_m$ 的简化剩余系．简化剩余系对模 $m$ 乘法封闭．简化剩余系中的每个元素都具有模 $m$ 意义下的唯一逆元，且其也在 ${\mathbb{Z}}_m^R$ 中．此二性质易由 $gcd$ 的线性组合意义证明．

简化剩余系的大小即 $\phi (m)$．这是 Euler 函数的其中一种定义．

设 $a\in {\mathbb{Z}}_m^R$，其在 ${\mathbb{Z}}_m$ 中存在逆元，故我们有 $$\mathrm{\forall }x,y\in {\mathbb{Z}}_m,x\equiv y(\mathrm{mod}m)⟺ax\equiv ay(\mathrm{mod}m)$$ 当限定 $x,y\in {\mathbb{Z}}_m^R$ 时，结合 ${\mathbb{Z}}_m^R$ 上述两个性质，就有 $${\mathbb{Z}}_m^R=a{\mathbb{Z}}_m^R=\{axmodm:x\in {\mathbb{Z}}_m^R\}$$ 故 $$a^{\phi (m)}\prod _{x\in {\mathbb{Z}}_m^R}x\equiv \prod _{x\in {\mathbb{Z}}_m^R}ax\equiv \prod _{x\in a{\mathbb{Z}}_m^R}x\equiv \prod _{x\in {\mathbb{Z}}_m^R}x⟹a^{\phi (m)}\equiv 1(\mathrm{mod}m)$$

$a^n$ 的取值 - $a$, $m$ 互质时

补充讨论 $a^{n}modm,n\in {\mathbb{Z}}_{+}$ 可取得哪些 ${\mathbb{Z}}_m$ 中的元素．

当 $gcd(a,m)=1$ 时，用 $gcd$ 的线性组合含义，有 $$\begin{array}{rl}gcd(a,m)=1& ⟹\mathrm{\exists }{x}_0,a{x}_0=1(\mathrm{mod}m)\\ & ⟹\mathrm{\exists }{x}_0,a^n{x}_0^n=1(\mathrm{mod}m)\\ & ⟹gcd(a^n,m)=1\end{array}$$ 故此时 $a^{n}modm$ 只可能取得 ${\mathbb{Z}}_m$ 的简化剩余系 ${\mathbb{Z}}_m^R$ 中的元素．

$a^n$ 的取值 - 一般情况⁹

当 $d_1=gcd(a,m)\ne 1$ 时，设 $a=d_1{a}_1,m=d_1{m}_1$，当 $n\ge 1$ 时， $$\begin{array}{rl}{a}^{n}modm& =d_1{a}_1{a}^{n-1}mod{d}_1{m}_1=d_1(a_1{a}^{n-1}mod{m}_1)\end{array}$$ 此时 $gcd(a_1,m_1)=1$，根据线性同余方程理论，括号内表达式的取值集合仅与 $a^{n-1}mod{m}_1$ 有关．此时若 $d_2=gcd(a,m_1)$ 仍不等于 $1$，则令 $a_1=d_2{a}_2,m_1=d_2{m}_2$．类似地继续化归，最终有 $d_{k+1}=gcd(a,m_k)=1$， $$a^{n}modm=d_1(a_1{d}_2(\dots a_{k-1}{d}_k(a_k{a}^{n-k}mod{m}_k)\dots )mod{m}_1)$$ 故当 $n\ge k$ 时，$a^{n}modm$ 的取值集合为 $$\{d_1{d}_2\dots d_k(a^{t}mod{m}_k):t\in {\mathbb{Z}}_{+}\}\subset d_1{d}_2\dots d_k{\mathbb{Z}}_{m_k}^R$$ 事实上，$\prod _{i=1}^k{d}_i$ 恰为 $gcd(a,m)$ 的各素因子在 $m$ 中全部成分的乘积．

原根¹⁰

需要注意的是，Euler 定理只给出了 ${\mathbb{Z}}_m$ 上元素阶的一个上界．换言之，其说明 $$\mathrm{\forall }a\in {\mathbb{Z}}_m^R⟹{\mathrm{ord}}_m(a)\mid \phi (m)$$ 而对于那些使得 $gcd(a,m)\ne 1$ 的元素 $a$，用 $gcd$ 的线性组合性质容易证明其永远不可能通过自乘变为 $1$．我们设定这类元素的阶为 $0$．

称 ${\mathbb{Z}}_m$ 上那些确能达到阶数上界 $\phi (m)$ 的元素是 ${\mathbb{Z}}_m$ 上的原根，用符号 $g$ 表示．

求阶和原根

• 求阶，用定理 ${\mathrm{ord}}_m(a)\mid \phi (m)$，求因子 + 快速幂即可 $O(\sqrt{m}\log m)$．
• 找原根，从小到大用上述求阶方法暴力即可．最小原根一般不会太大．
• 更快的方法请参考洛谷求原根模板题题解．

下面再介绍原根的两个定理．

必要性显然．充分性，反证出所有 $\phi (m)$ 的非平凡因子均不是 $g$ 的阶即可．

${\mathbb{Z}}_m$ 上的本原单位根构造

设 $a$ 是 ${\mathbb{Z}}_m$ 上一元素，当 $n\mid {\mathrm{ord}}_m(a)$，令 ${\omega }_n=a^{\frac{{\mathrm{ord}}_m(a)}{n}}$，容易发现 ${\omega }_n$ 就是 ${\mathbb{Z}}_m$ 上的 $n$ 次本原单位根．

使用数论中 Bézout 定理，${\omega }_n^k$ 是 ${\mathbb{Z}}_m$ 上的 $\frac{n}{gcd(n,k)}$ 次本原单位根．这是 ${\mathbb{Z}}_m$ 上本原单位根的消去引理．

遗憾的是，一般的 ${\mathbb{Z}}_m$ 中，本原单位根不一定满足求和引理¹¹ $\frac{1}{n}\sum _{k=0}^{n-1}{\omega }_n^{ik}=[n\mid i]$．这是因为，求和引理的证明中使用了等比数列求和和消去律，但一般的 ${\mathbb{Z}}_m$ 中消去律并不成立．

${\mathbb{Z}}_m$ 上的 NTT - $m$ 是质数

求和引理是保障 NTT 逆变换对应矩阵确为 $n^{-1}{\left(\begin{array}{c}{\omega }_n^{-ij}\end{array}\right)}_{(i,j)\in n\times n}$ 的关键性质．如果我们还想保留这一点，就必须要求 $m$ 是质数．

为保证 $n$ 在 ${\mathbb{Z}}_m$ 中存在逆元，必须要求 $gcd(n,m)=1$．已经要求 $m$ 是质数，故无需做出额外要求．

综上，逆变换矩阵确为 $n^{-1}{\left(\begin{array}{c}{\omega }_n^{-ij}\end{array}\right)}_{(i,j)\in n\times n}$ 的 $n$ 点 NTT 要求 $m$ 是质数，且 $n\mid {\mathrm{ord}}_m(a)$．已经要求 $m$ 是质数，故 ${\mathbb{Z}}_m$ 中必存在原根，不妨直接取 ${\mathbb{Z}}_m$ 的一个原根 $g$ 构造 $n$ 次本原单位根 ${\omega }_n=g^{\frac{\phi (m)}{n}}=g^{\frac{m-1}{n}}$，这样的构造支持至多 $m-1$ 次单位根的存在，显然是最优的选择．

${\mathbb{Z}}_m$ 上的 FFT 与卷积加速 - 998244353

为满足 FFT 对 $2^m$ 次本原单位根的需求，只需选择 $p=k{2}^m+1$ 型的奇质数 $p$，就可以在 ${\mathbb{Z}}_p$ 上支持规模至多为 $2^m$ 的 NTT/FFT 及其逆变换运行．

$998244353=2^{23}\times 7\times 17+1$，同时是一个质数，在 int 型中的单次加减操作不会溢出，是 OI/XCPC 计数题中不可多得的优秀模数¹²．

结合前述关于 ${\mathbb{Z}}_p[x]$ 上多点插值的讨论，通过 NTT/FFT 加速卷积运算的正确性得到完整证明．

${\mathbb{Z}}_m$ 上的 NTT - 对一般的 $m$

当 $m$ 不一定是质数时，NTT 逆变换不能表示为 $n^{-1}{\left(\begin{array}{c}{\omega }_n^{-ij}\end{array}\right)}_{(i,j)\in n\times n}$ 的形式．但这并不代表 NTT 不可逆．因此，尽管不再实用，一般 ${\mathbb{Z}}_m$ 上的 NTT 变换仍有讨论价值，即研究 ${\left(\begin{array}{c}{\omega }_n^{ij}\end{array}\right)}_{(i,j)\in n\times n}$ 的可逆性．

我们找到一篇有关该问题的参考文献[7, Secs. 3 and appendix B]，但尚不确定其证明的正确性．友情提示读者：${\mathbb{Z}}_m$ 不是域，甚至不是整环，故线性空间中矩阵的性质不能直接应用于 ${\mathbb{Z}}_m$ 上的矩阵，讨论可逆性时还需小心谨慎．

从分析到代数

后续数学内容导读 NTT 原理虽已非常”数学”，但也只是 DFT 在有限域上的一个实例．本节往后，我们要尝试为多项式系数位于复数域 $\mathbb{C}$ 上的 DFT、位于模 $p$ 剩余类域 ${\mathbb{Z}}_p$ 上的 NTT 及它们的快速算法 FFT 建立一个统一的数学框架．这需要我们剖析求值插值的基本原理，提炼出 FFT 算法成立的根本要求．

后续数学内容不再要求掌握．望同学们在纷繁的定理定义中抓住要旨，窥见抽象数学背后蕴藏的规律．熟悉高等代数和抽象代数的同学或会对某些内容感到熟悉．抽象地讨论 FFT 的资料并不多见，后续内容多为主讲人的新进探索，或有谬误，望不吝指正．

求值与插值 - 导言

本节将重新审视已经熟知的多项式，把抽象的、代数的多项式和具体的、分析的多项式函数区分开来．我们指出，多项式和多项式函数不同但关联紧密，形式幂级数与幂级数亦有此类联系．这些抽象的讨论将帮助我们剖析多项式求值插值的基本原理．

除常见代数书目（如[8]），也推荐参考 OI Wiki 的多项式基础[9]和 Wikipedia 的形式幂级数[10]．

多项式

设（无穷）序列 $\{a_n\}$ 是一个只有有限个非零项的序列，其元素均在环 $R$ 上．环 $R$ 上的多项式环[^14] $R[x]$ 是所有满足上述条件的序列构成的集合以及在它们之间定义的两种运算 $+,\times$ 的合称．此语境下，我们也将构成 $R[x]$ 的序列称为 $R$ 上的多项式．序列中的元素被称为多项式的系数．序列间定义的加法和乘法是普通序列升级为多项式的关键．

多项式环上的加法、乘法的定义已经为大家所熟知．系数所处的环保证了多项式加法和乘法的良定义，而在这两种运算下，$R$ 上多项式的集合也构成一个环的结构．

多项式

习惯上也会将多项式 $A$ 写为 $A(x)=\sum _{k=0}^{n-1}{a}_k{x}^k$ 的形式．其中未定元 $x$ 只是一个符号，仅为方便加、乘法的理解而使用，并无任何含义．

多项式 $A$ 的次数 $\mathrm{deg}A$ 定义为其最高非零项所处的位置下标．特别的，定义零多项式的次数为 $-\mathrm{\infty }$．对任意多项式 $A,B\in R[x]$，$\mathrm{deg}(A+B)\le max\{\mathrm{deg}A,\mathrm{deg}B\}$．当 $R$ 是整环¹⁴时，$A,B$ 的首项乘积非零，故 $\mathrm{deg}(A\ast B)=\mathrm{deg}A+\mathrm{deg}B$．

下面额外为多项式定义一种新的运算．多项式 $A$ 和多项式 $B$ 的复合 $A\circ B$ 定义为 $$A\circ B:=\sum _{k=0}^{+\mathrm{\infty }}{a}_k{B}^k$$ 其中 $B^k$ 代表 $k$ 个多项式 $B$ 的乘积．习惯上也会将多项式 $A$ 和多项式 $B$ 的复合写为 $A(B(x))=\sum _{k=0}^{+\mathrm{\infty }}{a}_k{B}^k(x)$ 的形式．由于多项式的次数有限，复合运算中仅包含了有限次多项式加法和乘法．

形式幂级数

形式幂级数定义与多项式的唯一区别是其不要求 $\{a_n\}$ 只有有限项非零．类似的，$R$ 上的形式幂级数环记为 $R[[x]]$．

由于涉及无限次运算，形式幂级数的复合运算需考虑环 $R$ 上的收敛问题．为回避此问题，一般规定复合右侧函数的常数项须为零．

带余除法

整环是无零因子的交换幺环．所谓无零因子，即环中任意元素 $a,b$ 满足 $ab=0⟹a=0\vee b=0$．所有域都是整环．当 $R$ 是整环时，$R[x]$ 也是整环¹⁵．整环上可以定义整除相关理论．

定义有带余除法的环被称为 Euclid 整环．域上的多项式环都是 Euclid 整环[8, 第 7 章第 2 节定理 3, p. 11]．值得注意的是，$R$ 的 Euclid 性不可传递至 $R[x]$¹⁶．

如在带余除法中保证除式是首一多项式¹⁷，则带余除法的良定义和进行过程也均可在整环上实现．

多项式函数 / 幂级数

刚刚强调，多项式 / 形式幂级数只是定义了加法和乘法的序列．现在介绍多项式函数和幂级数．它们不是序列，而是映射¹⁸．

多项式 $A$ 对应的环 $R$ 上的多项式函数 $A(x)$ 定义为映射¹⁹ $$A:R\to R,x\mapsto \sum _{k=0}^{+\mathrm{\infty }}{a}_k{x}^k$$ 其中 $\{a_n\}$ 只有有限项非零．环 $R$ 上的幂级数则无此限制²⁰．

与多项式 / 形式幂级数不同，这里的 $x\in R$ 代表函数的自变量，是会按环 $R$ 上运算法则参与运算的有意义的变量．

多项式函数的加法和乘法定义为函数的加法和乘法，即 $$\begin{array}{rl}A+B:x& \mapsto A(x)+B(x)\\ A\otimes B:x& \mapsto A(x)B(x)\end{array}$$

求值

多项式和多项式函数似乎在许多情况下有着平行的关系．下面介绍一个较直观的结论．

若环 $R$ 是交换环，则环 $R$ 上每个多项式都唯一确定一个 $R$ 上的多项式函数，且该映射保持加法和乘法运算²¹．形式化地可以记为 $$\begin{array}{rl}(A+B)(x)& =A(x)+B(x)\\ (A\ast B)(x)& =A(x)B(x)\end{array}$$ 可用于证明上述结论的定理证明见[8, 第 7 章第 1 节定理 4, p. 7]²²．由于多项式的复合由多项式加法和乘法组成，而多项式函数的复合（定义为函数的复合）亦可在环上平行地展开，故上述映射也保持了复合运算．

多项式和多项式函数的这一关系为多项式在任意点的求值操作提供了理论基础．在多项式函数的 $t$ 处求值也被称为用 $t$ 代入多项式 $A$，记为 $A(t)$．

插值

需要注意的是，前述结论的逆命题不一定成立，即环 $R$ 上的一个多项式函数不一定与 $R$ 上多项式一一对应．一个经典反例是有限域上多项式环 ${\mathbb{Z}}_p[x]$ 上的 $(x+1{)}^p\equiv x^p+1(\mathrm{mod}p)$．我们指出其成立的一个充分条件是 $R$ 是无限整环²³．

该逆命题的本质是通过多项式函数的所有函数值反过来确定多项式（系数）的过程．而如果在确定时只使用一部分函数值，就是所谓的多点插值过程．很多时候待求多项式的次数是已知的，这在相当程度上缩小了待定多项式的范围．我们指出，只要 $R$ 是整环，就可以通过多项式任意 $n$ 个不同位置的点值确定一个至多 $n-1$ 次的多项式．

证明的关键是用带余除法²⁴讨论多项式根与其一次因式的关系，即多项式余式定理或小 Bézout 定理[8, 第 7 章第 6 节定理 6, p. 35]．

需要强调，上述结论只能说明，只有那些确实可通过多项式 $n$ 点求值得到的点值组，才能与多项式建立保持加法和乘法运算的双射．对于任取的 $R$ 中的一个 $n$ 元组，则无法确定是否存在与之对应的多项式．

求值与插值的线性表示

到这里，多项式的求值与插值的线性表示已经呼之欲出了．若将交换环 $R$ 上的至多 $n-1$ 次的多项式看做模²⁵ $R^n$ 中的一个向量，则多项式多点求值的过程就是 $R^n$ 上的一个线性变换（记为 $V$），其矩阵表示正是 Vandermonde 矩阵．由先前对求值的讨论，我们还知道其保持模 $R^n$ 上的加法运算，并将结果不超过 $n-1$ 次的卷积运算转化为 $R^n$ 上的逐项乘法．

若将对 $R$ 的要求加强到域，注意到 Vandermonde 行列式在参数互不相同时非零，故此时 $V$ 可逆．事实上，Vandermonde 矩阵求逆的过程就是 Lagrange 插值，从而给出了将任意 $R^n$ 中向量还原为其对应多项式的构造方法²⁶．这与数论中中国剩余定理的构造有异曲同工之妙[11, section: A perspective from linear algebra]．

求值与插值 - 小结

本节的核心是多项式和多项式函数的区别与联系，两个方向的”确定”分别给出了多项式多点求值和多点插值的理论基础．

由于多项式和多项式函数这种若即若离的关系，往往在记号上也有意无意地混淆了它们，某些情况下加大了区分的难度．本篇使用的记号体系将尽量用单个字母 $A$ 表示多项式，用 $A(x)$ 表示多项式 $A$ 对应的多项式函数，用 $A(t)$ 表示用 $t$ 代入多项式 $A$，同时在需要区分的符号前加上适当的名词．

环上的 DFT - 导言

前面讨论了在多项式任意点处求值插值的基本原理，但 DFT/FFT 的运行只需在单位根处求值和插值．本节将进一步放宽对多项式环的限制，介绍定义在有主要单位根的环上的一般的 DFT 及其快速算法 FFT．

本节内容主要参考了 [12, Secs. 2, pp. 983–984] 和 [13]．

环上的单位根 - 两个定义

定义环 $R$ 上的 $n$ 次主要单位根（principal n-th root of unity）${\omega }_n$ 是满足如下条件的环 $R$ 中的一个元素： $${\omega }_n^n=1;\sum _{k=0}^{n-1}{\omega }_n^{ik}=0,i=1,2,\dots ,n-1$$ 由定义立得 $\sum _{k=0}^{n-1}{\omega }_n^{ik}=[n\mid i]n$，可见主要单位根的定义即求和引理本身．

定义环 $R$ 上的 $n$ 次本原单位根（primitive n-th root of unity）${\omega }_n$ 是满足如下条件的环 $R$ 中的一个元素： $${\omega }_n^n=1;{\omega }_n^k\ne 1,k=1,2,\dots ,n-1$$

本原单位根在许多情况下与主要单位根等价，但亦非完全相同．

环上的单位根 - 区别与联系

环上的单位根 - 区别与联系

环上的单位根 - 其它性质

分别根据定义和数论中的 Bézout 定理，容易证明主要单位根和本原单位根的消去引理：设 ${\omega }_n$ 是环上的 $n$ 次主要（本原）单位根，则 ${\omega }_n^k$ 是环上的 $\frac{n}{gcd(n,k)}$ 次主要（本原）单位根．

注意到 ${\left({\omega }_{2n}^n\right)}^2=1⟹({\omega }_{2n}^n+1)({\omega }_{2n}^n-1)=0$，故整环上的 $2n$ 次主要单位根或本原单位根 ${\omega }_{2n}$ 均满足 ${\omega }_{2n}^n=-1$．

环上的 $n$ 次主要单位根或本原单位根 ${\omega }_n$ 都存在逆元 ${\omega }_n^{-1}={\omega }_n^{n-1}$．

环上的 DFT

设 ${\omega }_n$ 是环 $R$ 上的一个 $n$ 次主要单位根，其对应的 $R^n$ 上的 $n$ 点 DFT 定义为线性映射（或矩阵） $$F={\left({\omega }_n^{ij}\right)}_{(i,j)\in n\times n}=\left(\begin{array}{cccc}1& 1& \dots & 1\\ 1& {\omega }_n& \dots & {\omega }_n^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{n-1}& \dots & {\omega }_n^{(n-1)(n-1)}\end{array}\right)$$ 若 $n1$ 在环 $R$ 中存在乘法逆元，则 $R^n$ 上的 DFT 可逆，其逆映射为²⁸ $$F^{-1}=(n1{)}^{-1}{\left({\omega }_n^{-ij}\right)}_{(i,j)\in n\times n}=(n1{)}^{-1}\left(\begin{array}{cccc}1& 1& \dots & 1\\ 1& {\omega }_n^{-1}& \dots & {\omega }_n^{-(n-1)}\\ ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{-(n-1)}& \dots & {\omega }_n^{-(n-1)(n-1)}\end{array}\right)$$

环上的 FFT

除无法在任意环上使用 ${\omega }_{2n}^n=-1$，现有的抽象已足够让我们写出与之前类似的推导 FFT 的过程．

设 $R$ 是任意环，${\omega }_{2n}$ 是环 $R$ 上的一个 $2n$ 次主要单位根．由主要单位根的消去引理，${\omega }_{2n}^2$ 是环 $R$ 上的 $\frac{2n}{gcd(2n,2)}=n$ 次主要单位根．方便起见，记 ${\omega }_n={\omega }_{2n}^2$．

考虑对模 $R^{2n}$ 中向量 $\mathit{a}={(a_0,a_1,\dots ,a_{2n-1})}^T$ 做 ${\omega }_{2n}$ 对应的 $2n$ 点 DFT 变换，得到向量 $\hat{\mathit{a}}={({\hat{a}}_0,{\hat{a}}_1,\dots ,{\hat{a}}_{2n-1})}^T$，其中 ${\hat{a}}_i=\sum _{k=0}^{2n-1}{a}_k{\omega }_{2n}^{ik}$．将其奇偶分项为 ${\hat{a}}_i=A_0(i)+A_1(i){\omega }_{2n}^i$，其中 $$\begin{array}{rl}{A}_0(i)& :=\sum _{k=0}^{n-1}{a}_{2k}{\omega }_n^{ik}=\sum _{k=0}^{n-1}{a}_{2k}{\omega }_{2n}^{2ik}\\ A_1(i)& :=\sum _{k=0}^{n-1}{a}_{2k+1}{\omega }_n^{ik}=\sum _{k=0}^{n-1}{a}_{2k+1}{\omega }_{2n}^{2ik}\end{array}$$