写在前面

1. 由于上一篇总结的版面限制，特开此文来记录 \(\text{OI}\) 中多项式类数学相关的问题。
2. 该文启发于Miskcoo的博客，甚至一些地方直接引用，在此特别说明；若文章中出现错误，烦请告知。
3. 感谢你的造访。

前置技能

多项式相关

形同 \(P(X)=a_0+a_1X+a_2X^2+\cdots+a_nX^n\) 的形式幂级数 \(P(X)\) 称为多项式。其中 \(\{a_i|i\in[0,n]\}\) 为多项式的系数； \(n\) 表示多项式的次数。

多项式的系数表示

对于 \(n\) 次多项式 \(P(X)\) 的系数 \(\{a_i|i\in[0,n]\}\) ，我们记向量 \(\vec a=(a_0,a_1,\cdots,a_n)\) 为 \(P(x)\) 的系数表示。

多项式的点值表示

对于一个 \(n\) 次多项式 \(P(X)\) ，我们由函数和方程的思想，对于函数 \(P(X)\) 的图像，我们只要确定了该图像上的 \(n+1\) 个点，那么 \(P(X)\) 是确定且唯一的。我们记这 \(n+1\) 个点组成的集合为 \(\{\left(x_i, P(x_i)\right)| i\in[0,n]\}\) 。那么该集合为多项式 \(P(x)\) 的点值表示。

一个多项式有不同的点值表示，但任意一个点值都能找出其唯一对应的多项式。

多项式的系数表示和点值表示能够互相转换。

复数相关

记 \(i=\sqrt{-1}\) 。我们把形如 \(a+bi\) （ \(a,b\) 均为实数）的数称为复数。其中 \(a\) 为虚数的实部， \(b\) 为虚数的虚部。

复数的意义

对于任意一个复数 \(a+bi\) ，都能用复平面上的唯一一个向量 \((a,b)\) 来表示。

复平面：即复数平面，由实轴作为 \(x\) 轴，虚轴作为 \(y\) 轴构成。

下文中出现的向量 \((a,b)\) 来代指虚数 \(a+bi\) 。

共轭复数：对于复数 \(z=a+bi\) ，对于另外一个复数 \(z'=a-bi\) ，我们称 \(z'\) 为 \(z\) 的共轭复数，记做 \(\overline{z}\) 。容易发现，一个复数与其共轭复数的实部相等，虚部互为相反数。

复数的辐角：我们可以将复数 \(z\) 写成 \(z=r\times(\cos\theta+i\sin\theta)\) ，其中 \(r\) 为复数 \(z\) 的模长， \(\theta\) 为复数 \(z\) 的辐角。一个复数有多个辐角，这些值相差 \(2\pi\) 。我们将 \(\theta\in[-\pi,\pi)\) 的 \(\theta\) 叫做辐角的主值。指数形式： \(z=r\times(\cos\theta+i\sin\theta)=re^{i\theta}\) 。

复数的基本运算

对于复数 \((a,b)\) 和 \((c,d)\) 。

满足加法法则 \((a,b)\pm(c,d)=(a\pm c,c\pm d)\) 。

满足乘法法则 \((a,b)\times(c,d)=(ac-bd,bc+ad)\) 。

对于除法，只需分子分母同乘分母的共轭复数，将分母实数化，分子做复数乘法即可。

复数乘法的几何意义：模长相乘，幅角相加。

证明：
对于两个复数 \(z_1=r_1\times(\cos\theta_1+i\sin\theta_1),z_2=r_2\times(\cos\theta_2+i\sin\theta_2)\) 。
\[\begin{aligned}z_1\times z_2&=r_1\times(\cos\theta_1+i\sin\theta_1)\times r_2\times(\cos\theta_2+i\sin\theta_2)\\&=r_1r_2\times((\cos\theta_1\cos\theta_2-\sin\theta_1\sin\theta_2)+i(\sin\theta_1\cos\theta_2+\cos\theta_1\sin\theta_2))\\&=r_1r_2\times(\cos(\theta_1+\theta_2)+i\sin(\theta_1+\theta_2))\end{aligned}\]
得证。

单位根

\(n\) 次单位根是指能够满足方程 \(z^n=1\) 的复数，这些复数一共有 \(n\) 个它们都分布在复平面的单位圆上，并且构成一个正 \(n\) 边形，它们把单位圆等分成 \(n\) 个部分。

根据复数乘法相当于模长相乘，幅角相加就可以知道， \(n\) 次单位根的模长一定是 \(1\) ，幅角的 \(n\) 倍是 \(0\) 。

这样， \(n\) 次单位根也就是

\[e^{\frac{2\pi ki}{n}}, k = 0, 1, 2, \cdots, n - 1\]

再根据欧拉公式

\[e^{\theta i}=\cos\theta + i\sin\theta\]

就可以知道 \(n\) 次单位根的算术表示

如果记 \(\omega_n=e^{\frac{2\pi i}{n}}\) ，那么 \(n\) 次单位根就是 \(\omega_n^0, \omega_n^1, \cdots, \omega_n^{n-1}\) 。

代码相关

给出此文出现的代码中的一些宏定义

#define dob complex<double>

const double pi = acos(-1.0);

const int mod = 998244353;

多项式乘法

给定两个多项式 \(A(x),B(x)\)

\[A(x) = \sum_{i=0}^na_ix^i = a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0 \\ B(x) = \sum_{i=0}^nb_ix^i = b_nx^n+b_{n-1}x^{n-1}+\cdots+b_1x+b_0\]

将这两个多项式相乘得到 \(C(x) = \sum_{i=0}^{2n}c_ix^i\) ，在这里

\[c_i=\sum_{j+k=i,0\leq j,k\leq n}a_jb_k\]

如果一个个去算 \(c_i\) 的话，要花费 \(O(n^2)\) 的时间才可以完成，但是，这是在系数表示下计算的，如果转换成点值表示，知道了 \(A(x),B(x)\) 的点值表示后，由于点数是 \(O(n)\) ，就可以直接将其相乘，在 \(O(n)\) 的时间内得到 \(C(x)\) 的点值表示。

由于 \(C(x)\) 的次数为 \(2n\) ，所以我们可以在 \(A(x),B(x)\) 上取 \(2n+1\) 个点，便于唯一确定 \(C(x)\) 。

如果能够找到一种有效的方法帮助我们在多项式的点值表示和系数表示之间转换，我们就可以快速地计算多项式的乘法了，快速傅里叶变换就可以做到这一点。

快速傅里叶变换

由刚才的设想，我们只需按序进行下述三次操作，即可降低多项式乘法的复杂度：

1. 将 \(A(x),B(x)\) 的系数表达转为点值表达；
2. 将 \(A(x),B(x)\) 的点值表达相乘；
3. 将得到的结果转为系数表达。

我们已知 2. 是 \(O(n)\) 。剩下的 1.和3. 可以用快速傅里叶变换来实现 \(O(n\log_2 n)\) 的变换。

DFT

\(\text{DFT}\) 是实现上述的 1. 过程的。它是一个基于分治策略的算法。

具体的思想是将多项式的 \(n\) 个系数通过变换变成 \(n\) 个点值。

对于一个多项式 \(A(x)=\sum_{i=0}^{n-1}a_ix^i\) ，我们先将 \(n\times 2\) ，理由上面说了，因为为了确定 \(C(x)\) 需要多取一些点。为了方便之后的处理，我们继续将 \(n\) 增大使得 \(n=2^m\) ，找到这个最小的 \(n\) ，并将不足的项的系数变为 \(0\) 。

接着将 \(n\) 个 \(n\) 次单位根 \(\omega_n^0, \omega_n^1, \cdots, \omega_n^{n-1}\) 代入 \(A(x)\)

\[A(\omega_n^k) = \sum_{i=0}^{n-1}a_i\omega_n^{ki} , k\in[0,n)\]

我们发现，这样计算出多项式 \(A(x)\) 的点值表示 \(\{\left(\omega_n^k, A(\omega_n^k)\right)| k\in[0,n)\}\) 的复杂度仍是 \(O(n^2)\) 的，如何将这个复杂度优化就是整个算法的关键。具体优化，就要利用单位根的性质了。

第一步，是将点集中每一项按照指数的奇偶分类：

\[\begin{eqnarray*} A(\omega_n^k) &=& \sum_{i=0}^{n-1}a_i\omega_n^{ki} \\ &=& \sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_n^{2ki}+\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_n^{2ki} \\ \end{eqnarray*}\]

但是这样分类后没什么用，因为对于每个系数 \(a_i\) ，它被计算的次数依旧是 \(n\) 次，因为对于每个 \(k\) 都要与 \(a_i\) 乘一次。

我们在回到上面这个式子，试着将它变形

\[A(\omega_n^k)=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\left(\omega_n^{ki}\right)^2+\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\left(\omega_n^{ki}\right)^2\]

注意到的是

\[\omega_n^2=\left(e^{\frac{2\pi i}{n}}\right)^2=e^{\frac{2\pi i}{\frac{n}{2}}}=\omega_{\frac{n}{2}}\]

这个等式可以直接由上述式子推出，但是我们可以去直观感受一下为什么会这样。
对于一个复数 \(z\) ，它的平方 \(z^2\) 依旧满足复数乘法运算的规律：模长相乘，幅角相加。
那么对于一个单位根 \(\omega_n=e^{\frac{2\pi i}{n}}\) ，它的平方模长不变，辐角翻倍。
我们想到一个复数有多个辐角，并且差为 \(2\pi\) 。既然如此，那么对于 \(n\) 次单位根 \(\omega_n\) ，与另一个 \(n\) 次单位根 \(\omega_n^{\frac{n}{2}}\) 辐角差 \(\pi\) 。平方后辐角翻倍那么相差 \(2\pi\) ，所以说这两个 \(n\) 次单位根的平方是相同的。
同时我们可以得到
\[\omega_n^{\frac{n}{2}+k} = \omega_n^{\frac{n}{2}}\cdot \omega_n^k = -\omega_n^k\]
因为 \(\omega_n^{\frac{n}{2}}=e^{\pi i}=-1\) ，得证。

即原式变成了

\[A(\omega_n^k)=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}+\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki}\]

不过这时是 \(k<\frac{n}{2}\) 时才成立的。由之前得到的式子 \(\omega_n^{\frac{n}{2}+k} = \omega_n^{\frac{n}{2}}\cdot \omega_n^k = -\omega_n^k\) 其余部分则应该满足：

\[\begin{eqnarray*} A(\omega_n^{k+\frac{n}{2}}) &=& \sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}+\omega_n^{k+\frac{n}{2}}\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki} \\ &=&\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}-\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki} \end{eqnarray*}\]

这样对于每个 \(a_i\) ，代入的值变成了 \(1, \omega_{\frac{n}{2}}, \omega_{\frac{n}{2}}^2, \cdots, \omega_{\frac{n}{2}}^{\frac{n}{2}-1}\) ，问题变成了两个规模减半的子问题，只要递归下去计算就可以了，复杂度是 \(O(n\log_2 n)\) 。

IDFT

刚才说了， \(\text{DFT}\) 是将系数表示转为点值表示。而我们现在需要解决的是将一个多项式的点值表示转为系数表示。 \(\text{IDFT}\) 则可实现这一过程。 \(\text{IDFT}\) 是 \(\text{DFT}\) 的逆过程。

考虑到一个本质的问题，如何基础地将点值表示转为系数表示，显然就是解如下的一个方程组：

\[\begin{equation*} \left\{ \begin{array}{ccccccccc} a_0(\omega_n^0)^{0}&+&\cdots&+&a_{n-2}(\omega_n^0)^{n-2}&+&a_{n-1}(\omega_n^0)^{n-1}&=&A(\omega_n^0) \\ a_0(\omega_n^1)^{0}&+&\cdots&+&a_{n-2}(\omega_n^1)^{n-2}&+&a_{n-1}(\omega_n^1)^{n-1}&=&A(\omega_n^1) \\ \vdots & & \vdots & &\vdots& & \vdots & & \vdots\\ a_0(\omega_n^{n-1})^{0}&+&\cdots&+&a_{n-2}(\omega_n^{n-1})^{n-2}&+&a_{n-1}(\omega_n^{n-1})^{n-1}&=&A(\omega_n^{n-1}) \end{array} \right. \end{equation*}\]

写成矩乘形式就是：

\[\begin{equation} \begin{bmatrix} (\omega_n^0)^0 & (\omega_n^0)^1 & \cdots & (\omega_n^0)^{n-1} \\ (\omega_n^1)^0 & (\omega_n^1)^1 & \cdots & (\omega_n^1)^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{n-1})^0 & (\omega_n^{n-1})^1 & \cdots & (\omega_n^{n-1})^{n-1} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix} = \begin{bmatrix} A(\omega_n^0) \\ A(\omega_n^1) \\ \vdots \\ A(\omega_n^{n-1}) \end{bmatrix} \end{equation}\]

左上的第一个矩阵为 \(\mathbf V\) 。现在考虑构造下面这个矩阵 \(d_{ij}=\omega_n^{-ij}\)

\[\begin{equation*} \mathbf D = \begin{bmatrix} (\omega_n^{-0})^0 & (\omega_n^{-0})^1 & \cdots & (\omega_n^{-0})^{n-1} \\ (\omega_n^{-1})^0 & (\omega_n^{-1})^1 & \cdots & (\omega_n^{-1})^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 & \cdots & (\omega_n^{-(n-1)})^{n-1} \end{bmatrix} \end{equation*}\]

记 \(\mathbf E=\mathbf D \cdot \mathbf V\) ，容易发现

\[\begin{eqnarray*} e_{ij} &=& \sum_{k=0}^{n-1} d_{ik} v_{kj} \\ &=& \sum_{k=0}^{n-1} \omega_n^{-ik}\omega_n^{kj} \\ &=& \sum_{k=0}^{n-1} \omega_n^{k(j-i)} \end{eqnarray*}\]

而对于这个式子容易发现

1. 当 \(i=j\) 时， \(\omega_n^{k(j-i)}=\omega_n^{0}=1\) ，故 \(e_{ij}=\sum_{k=0}^{n-1}1=n\) 。
2. 当 \(i\neq j\) 时， \(e_{ij}=\sum_{k=0}^{n-1}\omega_n^{k(j-i)}=\frac{\omega_n^0(1-\omega_n^{n(j-i)})}{1-\omega_n^{j-i}}\) ，由于 \(\omega_n^n=1\) ，故 \(e_{ij}=0\) 。

所以单位矩阵 \(\mathbf I_n=\frac{1}{n}\mathbf E\) ，由于 \(\mathbf I_n=\mathbf V\cdot\mathbf V^{-1},\mathbf E=\mathbf D \cdot \mathbf V\) ， \(\mathbf V\cdot\mathbf V^{-1}=\frac{1}{n}\mathbf D \cdot \mathbf V\) ， \(\mathbf V^{-1}=\frac{1}{n}\mathbf D\) 。其中 \(\mathbf V^{-1}\) 为 \(\mathbf V\) 的逆矩阵。

在上述 \((1)\) 式中等式左右两边左乘 \(\frac{1}{n}\mathbf D\) 。能得到

\[\begin{equation*} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix} = \frac{1}{n} \begin{bmatrix} (\omega_n^{-0})^0 & (\omega_n^{-0})^1 & \cdots & (\omega_n^{-0})^{n-1} \\ (\omega_n^{-1})^0 & (\omega_n^{-1})^1 & \cdots & (\omega_n^{-1})^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 & \cdots & (\omega_n^{-(n-1)})^{n-1} \end{bmatrix} \begin{bmatrix} A(\omega_n^0) \\ A(\omega_n^1) \\ \vdots \\ A(\omega_n^{n-1}) \end{bmatrix} \end{equation*}\]

即 \(a_k=\frac{1}{n}\sum_{i=0}^{n-1}=\omega_n^{-ki}A(\omega_n^i)\) ，这样， \(\text{IDFT}\) 就相当于把 \(\text{DFT}\) 过程中的 \(\omega_n^i\) 换成 \(\omega_n^{-i}\) ，然后做一次 \(\text{DFT}\) ，之后结果除以 \(n\) 就可以了。

算法实现

递归实现

由上面的说明想必很容易模拟上述过程，来实现递归的 \(\text{DFT}\) 和 \(\text{IDFT}\) 。

为了方便阅读，我们在这重写上面的式子，对于 \(\text{DFT}\) ：

\[\begin{aligned}A(\omega_n^k)&=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}+\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki}\\ A(\omega_n^{k+\frac{n}{2}}) &=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}-\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki}\end{aligned}\]

void FFT(dob *A, int len, int o) {

// len 为当前递归区间长度； o 为识别因子，若 o=1 表示进行 DFT ，若为 -1 表示进行 IDFT

    if(len == 1) return;

    dob wn(cos(2.0*pi/len), sin(2.0*pi*o/len)), w(1,0), t;

//注意此处 wn 的初值

    dob A0[len>>1], A1[len>>1];

    for (int i = 0; i < (len>>1); i++) A0[i] = A[i<<1], A1[i] = A[i<<1|1];

    FFT(A0, len>>1, o); FFT(A1, len>>1, o);

    for(int i = 0; i < (len>>1); i++, w *= wn) {

        t = w*A1[i];

        A[i] = A0[i]+t;

        A[i+(len>>1)] = A0[i]-t;

    }

}

递归来实现 \(\text{FFT}\) 的显著的优点就是直观简洁，几乎就是按照式子模拟即可；不过缺点是常数巨大，在实战上毫不适用。

迭代实现

假设现在有 \(16\) 个数要进行 \(\text{DFT}\) 来看看递归的过程。

（图片转自Miskcoo）

在 \(\text{Step1} \rightarrow \text{Step2}\) 的过程中，按照奇偶分类，二进制位中最后一位相同的被分到同一组；

在 \(\text{Step2} \rightarrow \text{Step3}\) 的过程中，仍然按照奇偶，只不过不是按照数字的奇偶性，而是下标的奇偶性，二进制位中最后两位相同的才被分到同一组；

在 \(\text{Step3} \rightarrow \text{Step4}\) 的过程中，二进制位中最后三位相同的数字被分在同一组；

我们将所有数的二进制翻转，容易发现每次分在同一组内的都是前缀相同的，并且同一组内的数值是连续的。

由于迭代实现，我们可以先将原来的 \(A\) 数组的位置预处理成 \(B\) 数组——递归最下面一层（最后一步）的 \(A\) 的位置。

假设 reverse(i) 是将二进制位反转的操作，那么 \(A\) 和 \(B\) 之间有这样的关系 B[reverse(i)]=A[i] ，也就是说， B[i+1]=A[reverse(reverse(i)+1)] ， \(B\) 中第 \(i\) 项的下一项就是将 \(i\) 反转后加 \(1\) 再反转回来 \(A\) 中的那一项，所以现在要模拟的就是从高位开始的二进制加法。

我们可以处理一个数组 \(R\) ： R[i]=reverse(i) ，即 B[R[i]]=A[i] 。

对于 \(R\) 的预处理，可以 \(O(n)\) 递推出，其中 \(L\) 表示 \(n=2^L\) ：

for (int i = 0; i < n; i++) R[i] = (R[i>>1]>>1)|((i&1)<<(L-1));

显然满足 R[R[i]]=i ，所以对于一组 \(i,R_i\) ，只要交换一次位置即可。

既然已经处理最后的位置 \(B\) 。那么我们就可以枚举区间长度迭代计算了。

void FFT(dob *A, int o) {

    for (int i = 0; i < n; i++) if (i < R[i]) swap(A[i], A[R[i]]);

    for (int i = 1; i < n; i <<= 1) {

//枚举区间长度

        dob wn(cos(pi/i), sin(pi*o/i)), x, y;

        for(int j = 0; j < n; j += (i<<1)) {

//枚举区间左端点

            dob w(1, 0);

            for(int k = 0; k < i; k++, w *= wn) {

//枚举 k

                x = A[j+k]; y = w*A[j+i+k];

                A[j+k] = x+y;

                A[j+k+i] = x-y;

            }

        }

    }

}

快速数论变换

有些题目要求答案要对一个质数取模，用 \(\text{FFT}\) 难免会有缺陷了；毕竟在复数域计算难免有精度损失。

首先来看 \(\text{FFT}\) 中能在 \(O(n\log_2n)\) 时间内变换用到了单位根 \(\omega\) 的什么性质。

1. \(\omega_n^n=1\) ，有这样一个初值，方便后面的计算；
2. \(\omega_n^0, \omega_n^1, \cdots, \omega_n^{n-1}\) 是互不相同的，这样带入计算出来的点值才可以用来还原出系数；
3. \(\omega_n^2=\omega_{\frac{n}{2}}, \omega_n^{\frac{n}{2}+k}=-\omega_n^k\) ，这使得在按照指数奇偶分类之后能够把带入的值也减半使得问题规模能够减半。
4. \[\sum_{k=0}^{n-1} (\omega_n^{j-i})^k = \begin{eqnarray*} \left\{ \begin{aligned}0, ~~~&i \neq j\\ n, ~~~&i = j \end{aligned} \right. \end{eqnarray*}\]
   这样保证了能够使用相同的方法进行逆变换得到系数表示。

如果我们能够在数论域上找到这样的类似单位根的东西。就可以用相同的思想来进行简化运算。

这样我们引出了原根的概念。

原根

根据费马小定理我们知道，对于一个素数 \(p\) ，有下面这样的关系

\[a^{p-1} \equiv 1 \pmod p\]

这样类比单位根，可以满足“周期性”。

对于一个数 \(g\) 满足 \(g^0, g^1, \cdots, g^{p-2} \pmod p\) 互不相同，那么称 \(g\) 是 \(p\) 的原根。

令 \(n=2^k\) ，我们取素数 \(p = r\cdot 2^k + 1\) （满足这样形式的素数叫做费马素数，朴素的 \(\text{NTT}\) 是只在费马素数下适用的。），然后取 \(p\) 的原根 \(g\) ，然后我们令 \(g_n\equiv g^r\pmod p\) ，这样就能满足 \(g_n^0, g_n^1, \cdots, g_n^{n-1} \pmod p\) 互不相同，并且 \(g_n^n\equiv 1\pmod p\) 。

首先 \(g_n^n\equiv 1\pmod p\) 是显然的，因为 \(2^{rn}\geq r\cdot 2^n\) ，就会满足 \(r\cdot 2^n\mid 2^{r+n}\) ；
而对于 \(g_n^0, g_n^1, \cdots, g_n^{n-1} \pmod p\) 互不相同的证明，其实等价于证明 \(0r,1r,\cdots,(n-1)r \pmod{p-1}\) 互不相同，这个在之前的博客中有讲，就不赘述了。

于是这样就满足了上面所说的 1.2. 。

由于 \(p\) 是质数且 \(g_n^n \equiv 1 \pmod p\) ，所以 \(g_n^{\frac{n}{2}} \equiv \pm 1 \pmod p\) ，又由于 2. ，所以 \(g_n^{\frac{n}{2}} \equiv -1 \pmod p\) 。直接满足了性质 3. 。

对于 4. ，可以用相同的方法代入计算。不再赘述。

我们这样就找到了这样一系列满足条件的数 \(g_n\) 。

对于一个模数，需要去找到它的原根，似乎比较麻烦，但Miskcoo提供了一些常用的数值。查表就好了。

算法实现

实现和 \(\text{FFT}\) 类似，我们依旧用迭代实现。

void NTT(int *A, int o) {

    for (int i = 0; i < n; i++) if (i < R[i]) swap(A[i], A[R[i]]);

    for (int i = 1; i < n; i <<= 1) {

//枚举区间长度

        int gn = quick_pow(3, (mod-1)/(i<<1)), x, y;

        if (o == -1) gn = quick_pow(gn, mod-2);

        for (int j = 0; j < n; j += (i<<1)) {

//枚举区间左端点

            int g = 1;

            for (int k = 0; k < i; k++, g = 1ll*g*gn%mod) {

//枚举 k

                x = A[j+k], y = 1ll*g*A[j+k+i]%mod;

                A[j+k] = (x+y)%mod;

                A[j+k+i] = (x-y+mod)%mod;

            }

        }

    }

}

模数任意的解决方案

传说中的 \(\text{MTT}\) （快速毛爷爷变换），具体思想是取几个乘积大于 \(n(mod-1)^2\) 的费马素数作为模数，分别求出结果后用 \(crt\) 合并就好了。

应用

快速卷积

对于两个定义在 \(\mathbb{N}\) 上的函数 \(f(n),g(n)\) ，定义 \(f\) 和 \(g\) 卷积为 \(f\otimes g\)

\[(f \otimes g)(n) = \sum_{i=0}^n f(i)g(n-i)\]

容易发现两个数论函数的卷积其实就是两个多项式的乘积。

可以用 \(\text{FFT}\) 或 \(\text{NTT}\) 优化。

更一般地，对于 \(h(k) = \sum_{i=0}^n f(i)g(i+k)\) ，我们可以设 \(f'(x)=f(n-x)\) ，显然 \(h(k) = \sum_{i=0}^n f'(n-i)g(i+k)\) ，也是一个卷积的形式。

多项式求逆

基本概念

对于多项式 \(A(x),B(x)\) ，存在唯一的 \(Q(x),R(x)\) 满足 \(A(x)=B(x)Q(x)+R(x)\) ，其中 \(R\) 的次数小于 \(B\) 的次数，我们称 \(Q(x)\) 为 \(B(x)\) 除 \(A(x)\) 的商， \(R(x)\) 为 \(B(x)\) 除 \(A(x)\) 的余数，可以记作

\[A(x) \equiv R(x) \pmod {B(x)}\]

对于一个多项式 \(A(x)\) ，如果存在 \(B(x)\) 满足 \(B\) 的次数小于等于 \(A\) 的次数，并且

\[A(x)B(x) \equiv 1 \pmod {x^n}\]

那么称 \(B(x)\) 为 \(A(x)\) 在 \(\mod x^n\) 意义下的逆元，记作 \(A^{−1}(x)\) 。

对于多项式 \(A(x)\) ，对于 \(A(x)\mod x^n\) 直观的解释是提出 \(A(x)\) 中的 \(x^0\sim x^{n-1}\) 次项。

求解方法

考虑如何求 \(A^{-1}(x)\) 。

1. 当 \(n=1\) 时，容易发现 \(A(x)\equiv c\pmod x\) 此时 \(c\) 是一个常数，故 \(A^{-1}(x)\) 也是一个常数，即 \(c^{-1}\) 。
2. 当 \(n>1\) 时，记 \(B(x)=A^{-1}(x)\) ，此时应该满足
   \[\begin{equation}A(x)B(x)\equiv 1\pmod {x^n}\end{equation}\]
   假设在 \(\mod x^{\lceil \frac{n}{2} \rceil}\) 意义下 \(A(x)\) 的逆元是 \(B'(x)\) 并且我们已经求出，那么
   \[\begin{equation}A(x)B'(x) \equiv 1 \pmod {x^{\lceil \frac{n}{2} \rceil}} \end{equation}\]
   再将 \((2)\) 放在 \(\mod x^{\lceil \frac{n}{2} \rceil}\) 意义下
   \[\begin{equation}A(x)B(x)\equiv 1\pmod {x^{\lceil \frac{n}{2} \rceil}}\end{equation}\]
   由 \((3)-(4)\) 得到
   \[B(x) - B'(x) \equiv 0 \pmod {x^{\lceil \frac{n}{2} \rceil}}\]
   两边同时平方
   \[B^2(x) - 2B'(x)B(x) + B'^2(x) \equiv 0 \pmod {x^n}\]
   解释一下平方后为什么模的 \(x^{\lceil \frac{n}{2} \rceil}\) 也会平方。
   因为，左边多项式在 \(\mod x^n\) 意义下为 \(0\) ，那么就说明其 \(0\) 到 \(n-1\) 次项系数都为 \(0\) ，平方了之后，对于第 \(0\leq i\leq 2n-1\) 项，其系数 \(a_i\) 为 \(\sum_{j=0}^i a_ja_{i-j}\) ，很明显 \(j\) 和 \(i-j\) 之间必然有一个值小于 \(n\) ，因此 \(a_i\) 必然是 \(0\) ，也就是说平方后在 \(\mod x^{2n}\) 意义下仍然为 \(0\) 。
   这时我们只要在等式两边同乘上 \(A(x)\) ，移项得
   \[B(x) \equiv 2B'(x) - A(x)B'^2(x) \pmod {x^n}\]
   这样就可以得到 \(\mod x^n\) 意义下的逆元了，利用 \(\text{FFT}\) 加速之后可以做到在 \(O(n\log_2n)\) 时间内解决当前问题。

由主定理，最后总的时间复杂度也就是 \(O(n\log_2n)\) 的。

顺便一提，由这个过程可以看出，一个多项式有没有逆元完全取决于其常数项是否有逆元。

算法实现

#include <bits/stdc++.h>

using namespace std;

const int mod = 998244353;

const int N = 4e5;

int n, a[N+5], b[N+5], L, R[N+5], len, tmp[N+5];

int quick_pow(int a, int b) {

    int ans = 1;

    while (b) {

        if (b&1) ans = 1ll*ans*a%mod;

        b >>= 1, a = 1ll*a*a%mod;

    }

    return ans;

}

void NTT(int *A, int o) {

    for (int i = 0; i < len; i++) if (i < R[i]) swap(A[i], A[R[i]]);

    for (int i = 1; i < len; i <<= 1) {

        int gn = quick_pow(3, (mod-1)/(i<<1)), x, y;

        if (o == -1) gn = quick_pow(gn, mod-2);

        for (int j = 0; j < len; j += (i<<1)) {

            int g = 1;

            for (int k = 0; k < i; k++, g = 1ll*g*gn%mod) {

                x = A[j+k], y = 1ll*A[j+k+i]*g%mod;

                A[j+k] = (x+y)%mod;

                A[j+k+i] = (x-y+mod)%mod;

            }

        }

    }

}

void poly_inv(int *A, int *B, int deg) {

// deg 表示多项式的度

    if (deg == 1) {B[0] = quick_pow(A[0], mod-2); return; }

    poly_inv(A, B, (deg+1)>>1);

    for (L = 0, len = 1; len <= (deg<<1); len <<= 1) ++L;

// A*B 的度为 deg^2

    for (int i = 0; i < len; i++) R[i] = (R[i>>1]>>1)|((i&1)<<(L-1));

    for (int i = 0; i < deg; i++) tmp[i] = A[i];

    for (int i = deg; i < len; i++) tmp[i] = 0;

    for (int i = (deg+1)>>1; i < len; i++) B[i] = 0;

//注意将高次项系数补为 0

    NTT(tmp, 1), NTT(B, 1);

    for (int i = 0; i < len; i++)

        B[i] = 1ll*B[i]*((2ll-1ll*tmp[i]*B[i]%mod+mod)%mod)%mod;

    NTT(B, -1); int inv = quick_pow(len, mod-2);

    for (int i = 0; i < len; i++) B[i] = 1ll*B[i]*inv%mod;

}

void work() {

    scanf("%d", &n); for (int i = 0; i < n; i++) scanf("%d", &a[i]);

    poly_inv(a, b, n); for (int i = 0; i < n; i++) printf("%d ", b[i]);

}

int main() {work(); return 0; }

求第二类斯特林数

第二类斯特林数

定义

将 \(n\) 个有区别的球放入 \(m\) 个无区别的盒子中非空的方案数，记为 \(S(n,m)\) 或 \(\begin{Bmatrix}n\\m\end{Bmatrix}\) 。

递推式

1. \(\begin{Bmatrix}i\\0\end{Bmatrix}=0,i\in\mathbb{N_+}\) ， \(\begin{Bmatrix}0\\0\end{Bmatrix}=1\)
2. \(\begin{Bmatrix}n\\m\end{Bmatrix}=\begin{Bmatrix}n-1\\m-1\end{Bmatrix}+m\begin{Bmatrix}n-1\\m\end{Bmatrix}\)

① 边界情况显然；
② 考虑第 \(n\) 个球如何放：放在之前的盒子里面，则共 \(m\) 方法；或新开一个盒子。

通项公式

\[S(n,m)=\frac{1}{m!} \sum _{k=0}^m (-1)^k{m\choose k}(m-k)^n\]

证明的话大致就是容斥原理， \(k\) 枚举有多少个集合是空的，每种情况有 \(m\choose k\) 种空集情况，\(n\) 个元素可以放进非空的 \(m-k\) 个集合中。这样求出来的答案是有序的，所以我们除以 \(m!\) 使得其变为无序。

\(\text{NTT}\) 优化

由第二类斯特林数的通项公式，我们将组合数拆开，得到

\[S(n,m)=\sum _{k=0}^m \frac{(-1)^k}{k!}\frac{(m-k)^n}{(m-k)!}\]

记多项式

\[\begin{aligned}C(x)&=\sum_{i=0}^\infty S(n,i)x^i\\A(x)&=\sum_{i=0}^\infty \frac{(-1)^i}{i!}x^i\\B(x)&=\sum_{i=0}^\infty\frac{i^n}{i!}x^i\end{aligned}\]

那么 \(C(x)=A(x)B(x)\) 。 \(\text{NTT}\) 优化即可。

快速沃尔什变换

我们回到多项式乘法的概念：我们已知多项式 \(A(x),B(x)\) 。要求多项式 \(C(x)=A(x)B(x)\) 。其中

\[c_i=\sum_{j+k=i}a_jb_k\]

似乎求和式下面的 j+k=i 比较单调。我们试着将 + 换成其他的符号。

考虑如何求

\[c_i=\sum_{j\oplus k=i}a_jb_k\]

其中 \(\oplus\) 为按位运算符号。包括常用的 xor and or ，即“按位异或”，“按位与”，“按位或”。

由于这三种卷积具有相似性，这里仅举 xor 为例，其余两种可以类比推出。

\(xor\) 卷积

注意下文中的 \(\oplus\) 表示“按位异或”。

我们类比 \(\text{FFT}\) 的过程。对于 \(C(x)=A(x)B(x)\) ，那么满足

\[\text{DFT}(A(x))_i\times\text{DFT}(B(x))_i=\text{DFT}(C(x))_i\]

其中 \(\text{DFT}\) 表示多项式系数表示转点值表示的过程。

我们考虑是否能也构造一个变换 \(\text{tf}\) ，使得 \(C(x)=A(x)\oplus B(x)\) 满足

\[\text{tf}(A(x))_i\times\text{tf}(B(x))_i=\text{tf}(C(x))_i\]

由于是“按位异或”，我们不妨先考虑只含一位的情况，此时多项式只有两项分别为 \(0\) 和 \(1\) 。记 \(A_0=0,A_1=1\) 。 \(B,C\) 类似。那么

\[\begin{aligned}\text{tf}(A)&=<A_0+A_1,A_0-A_1>\\\text{tf}(B)&=<B_0+B_1,B_0-B_1>\\\text{tf}(C)&=<C_0+C_1,C_0-C_1>\end{aligned}\]

至于为什么是这种形式，可以结合式子

\[\text{tf}(A(x))_i\times\text{tf}(B(x))_i=\text{tf}(C(x))_i\]

得出 \[\begin{cases}C_0&=A_0B_0+A_1B_1\\C_1&=A_0B_1+A_1B_0\end{cases}\]

显然这是满足异或卷积的表达式的，成立。

推广到一般的情况，当 \(A\) 的长度为 \(2^k\) 时，我们记 \(A_0\) 为前 \(2^{k-1}\) 位， \(A_1\) 为后 \(2^{k-1}\) 位。用数学归纳法证明。容易发现 \(A_0\) 中每一项的最高位为 \(0\) ， \(A_1\) 中每一项的最高位为 \(1\) 。

只要证明满足

\[\begin{aligned}\text{tf}(A)&=<\text{tf}(A_0)+\text{tf}(A_1),\text{tf}(A_0)-\text{tf}(A_1)>\\\text{tf}(B)&=<\text{tf}(B_0)+\text{tf}(B_1),\text{tf}(B_0)-\text{tf}(B_1)>\\\text{tf}(C)&=<\text{tf}(C_0)+\text{tf}(C_1),\text{tf}(C_0)-\text{tf}(C_1)>\end{aligned}\]

时， \(\text{tf}(A(x))_i\times\text{tf}(B(x))_i=\text{tf}(C(x))_i\) 依旧成立即可。

我们依旧暴力拆开，得到

\[\begin{cases}\text{tf}(C_0)&=\text{tf}(A_0)\text{tf}(B_0)+\text{tf}(A_1)\text{tf}(B_1)\\\text{tf}(C_1)&=\text{tf}(A_0)\text{tf}(B_1)+\text{tf}(A_1)\text{tf}(B_0)\end{cases}\]

由于异或每一位是独立，而这里如果我们把 \(C\) 按照最高位为 \(0\) 或 \(1\) 分成两部分，最高位的异或和其它位不相关。而 \(\text{tf}\) 已经将除最高位外的其他所有位处理好了。显然是满足条件的。注意： \(\text{tf}\times\text{tf}\) 是按位乘的，而不是卷积。

似乎我们已经做好了类似“系数转点值”的过程；考虑逆过程 \(\text{utf}\) 。依旧用同样的方法验证，先考虑只含一位的情况。

\[\text{utf}(A)=\left<\frac{A_0+A_1}{2},\frac{A_0-A_1}{2}\right>\]

由于 \(\text{uft}\) 是多项式，注意是“卷积”的形式。

得到 \(\begin{cases}C_0&=A_0B_0\\C_1&=A_1B_1\end{cases}\) 。显然满足。

接着就直接数归来证即可。方法类似之前证明 \(\text{tf}\) 的过程。

结论（三种卷积求法）

正向 \(\text{tf}\)

1. \(xor\) 卷积： \(\text{tf}(A)=<A_0+A_1,A_0-A_1>\)
2. \(and\) 卷积： \(\text{tf}(A)=<A_0+A_1,A_1>\)
3. \(or\) 卷积： \(\text{tf}(A)=<A_0,A_0+A_1>\)

逆向 \(\text{utf}\)

1. \(xor\) 卷积： \(\text{utf}(A)=\left<\frac{A_0+A_1}{2},\frac{A_0-A_1}{2}\right>\)
2. \(and\) 卷积： \(\text{utf}(A)=<A_0-A_1,A_1>\)
3. \(or\) 卷积： \(\text{utf}(A)=<A_0,A_1-A_0>\)

算法实现

这里仅提供 \(xor\) 卷积的模板，其它情况类似。

值得注意的是 \(\text{FWT}\) 与 \(\text{NTT},\text{FFT}\) 并不完全相似。只是用了类比的思想得到 \(\text{FWT}\) 的 \(\text{tf}\) 和 \(\text{utf}\) 过程，本质上是不同的。

所以说代码实现也只是借用了迭代的思想等，一些操作如“交换初始位置”，“逆变换后除以 \(n\) ”，是不需要的。

void FWT(int *A, int o) {

    for (int i = 1; i < n; i <<= 1)

        for (int j = 0; j < n; j += (i<<1))

            for (int k = 0; k < i; k++) {

                int x = A[k+j], y = A[k+j+i];

                A[k+j] = (x+y)%mod, A[k+j+i] = (x-y+mod)%mod;

                if (o == -1)

                    A[k+j] = 1ll*A[k+j]*inv2%mod,

                    A[k+j+i] = 1ll*A[k+j+i]*inv2%mod;

            }

}