当前位置：首页 > news >正文

简单多项式

news 2025/10/24 6:14:56

前言：本文按照笔者自身理解撰写，力求通透易懂，如有不严谨处请多担待，有问题欢迎反馈。

复数：

我们都知道 $x ^ 2 = 1$ 的解是：$x_1 = 1,x_2 = -1$ 那么方程 $x ^ 4 = 1$ 的解是什么呢？

首先一个 $N$ 次方程如果有 $N$ 个解，所以这里应该有四个解。

因为 $x ^ 4 = (x ^ 2) ^ 2$ 所以 $x ^ 2 = 1$ 或 $x ^ 2 = -1$ 这时我们会发现 $x = \sqrt{-1}$ 我们称这样的 $x$ 叫做 $i$，同理 $i$ 也满足 $i ^ 2 = -1$。

这里引入复数：对于任意的数，都可以表示成 $a + bi$ 这里 $a,b\in R$ 并且实数是完全被复数包含在内的。

我们知道，对于任意实数 $x\in R$ 其在数轴上都能被唯一地表示出来。那么对于复数，我们应该用什么样的坐标系，才能让他们都能被唯一地表示出来呢？

考虑每个 $a + bi$ 我们将其常数项和 $i$ 的系数组成一个二元组，发现对于每种二元组都能唯一地对应一个复数，或者说，每个这样的二元组都和其对应的复数组成双射关系，那么我们就可以把它们放到二维坐标系中。

一些概念：

我们把复数 $a + bi$ 对应的坐标 $(a,b)$ 所在的坐标系中的横轴表示实轴，纵轴表示虚轴。
如果我们把原点 $O$ 和一个复数 $a + bi$ 对应的坐标 $(a,b)$ 连边，那么这条边的长度就称作这个复数的模长。
如果我们把原点 $O$ 和一个复数 $a + bi$ 对应的坐标 $(a,b)$ 连边，那么这条边和实轴的正方向所夹的角称作幅角。

由于每个复数所对应的二元组在二维坐标系中，我们可以将它和向量的知识联系在一起。

复数运算：

加法运算，设我们有两个复数 $a + bi$ 和 $c + di$ 那么运算法则就是对应位置（即：常数相加，$i$ 的系数相加）相加，也就是 $a + bi + c + di \rightarrow a + c + bi + di \rightarrow a + c + (b + d)i$
减法运算，其实就是加法运算的逆，可以把 $a + bi - (c + di)\rightarrow a + bi + (-c - di)$ 即 $a + bi - (c + di) \rightarrow a - c + (b - d)i$
乘法运算，其实就是 $(a + bi)\times (c + di) \rightarrow ac + bdi^2 + bci + adi$ 因为 $i ^ 2 = -1$，所以原式就是：$(a + bi)\times (c + di) \rightarrow ac - bd + (bc + ad)i$
除法运算，其实就是 $\frac{a + bi}{c + di}$ 可以分子分母同乘 $c - di$ 从而分子变成平方差，即 $\frac{a + bi}{c + di} \times \frac{c - di}{c - di} \rightarrow \frac{ac - adi + bci -bdi^2}{c ^ 2 - d^2i^2}\rightarrow \frac{ac + bd + (bc - ad) i}{c ^ 2 + d ^ 2}$
复数的共轭：实部相同，虚部相反(根据实轴对称)。
复数相乘的时候模长相乘，幅角相加。

单位圆：

和三角函数一样，就是以 $O$ 为圆心，半径为一个单位长度的圆。
比如这样：

不难发现，单位圆上的每一个点模长都是 $1$
为了方便书写，我们下文称 $|x|$ 为 $x$ 的模长。

快速傅里叶变换思想：

引理：一个 $N$ 次多项式可以由 $N + 1$ 个 $x$ 互不相同的点值唯一确定出来，比如一次函数需要两个点，二次函数需要三个点。

一些符号的说明：设当前有两个多项式 $f,g$，我们称 $fg$ 表示多项式 $f$ 和 $g$ 相乘，即，若 $h = fg$ 则 $h_i = \sum_{j + k = i} f_j\times g_k$，定义 $deg_f$ 表示多项式 $f$ 次数。

首先暴力做法就是对于 $0\leq i\leq deg_f,0\leq j\leq deg_g, h_{i + j}\rightarrow h_{i + j} + f_i\times g_j$
时间复杂度是 $O(n^2)$ 的。

我们设多项式 $f(x) = a_0 + a_1x + a_2x^2 + ... a_nx^n$ 中系数为 $<a_0,a_1,....a_n>$ 我们发现直接对系数进行合并是难以优化的，但是如果转成点值呢？假设我们已经把 $f,g$ 分别转成点值序列 $A$ 和 $B$ 那么点值相乘是 $O(n)$ 的，即 $C = A\times B,\forall i\in [0,\max(deg_f,deg_g)],C_i = A_i\times B_i$ 最后我们只需要再进行一次多项式转点值的逆，也就是点值转多项式就好了，所以我们把重点放到了如何知道一个多项式 $f$ 求出 $deg_f + 1$ 个点的点值序列。

离散傅里叶变换 DFT

既然我们打算要把一个多项式 $f$ 转成一个点值序列 $A$ 那么我们就需要任意找 $deg_f + deg_g + 1$ 个点分别带入两个多项式，因为 $|fg| = |f| + |g| - 1$ 而 $|f| = deg_f + 1$ 其中 $|f|$ 表示多项式 $f$ 有多少个幂， $0$ 次幂也算幂。

但是想法是好的，显示是很残酷的，因为多项式多点求值需要使用 $FFT$ 实现，但是我们研究的是 $FFT$ 肯定是不能这样的，因此我们应该找一些特殊点。

伟大的傅里叶告诉我们，要从单位根上来考虑这个问题。

具体的，我们把单位圆均分成 $N$ 份，比如 $N = 8$ 的时候图长这样：

我们发现，这 $N$ 个根也是 $x ^ N = 1$ 的解，后文称作单位根，我们设实轴上的边为 $\omega_{n}^0$ 然后逆时针的过程中第一个遇到的红线就是 $\omega_{n}^1$ 同理，第 $i$ 个遇到的红线就是 $\omega_{n}^i$ 容易发现，由于我们有 $n$ 个单位根，所以也就是有 $\omega_{n}^0 .... \omega_{n}^{n - 1}$ 下面我们重点来讨论一下这些单位根的性质。

单位根性质：

$\omega_n^{0,1.....n - 1}$ 并不是几个独立的点，而是代表一个集合，因为如果在某一个点，它逆时针再走 $2\pi$ 步仍然可以到这个点，由此也能推出单位圆中的周期性 $T = 2\pi$ 或者说 $w_n^i = w_n^{i\bmod n},i\in Z$
$\omega_n^j \times \omega_n^i = \omega_n^{i + j}$ 这个就是刚刚复数运算中的复数相乘，幅角相加。
$(\omega_n^i)^j = \omega_n^{ij}$ 这个其实就是 $j$ 个 $w_n^i$ 相乘和 $2$ 一样的分析。
$\omega_{2n}^{2i} = \omega_n^i$ 类比分蛋糕的时候，就是原来是分了 $2n$ 块然后选了 $2i$ 块，现在是分了 $n$ 块，选了 $i$ 块，不难发现是一样的，注意：这里的“分”指的是均分。
如果 $n$ 是偶数，那么 $\omega_n^i = -\omega_n^{i + \frac{n}{2}}$ 因为 $n$ 是偶数，所以 $w_n^{\frac{n}{2}}$ 就相当于逆时针走了 $180$ 度，自然就是变成了原来中心对称的位置了。

准备好了吗，即将开始最重要的部分！

考虑当前有个 $n$ 次多项式 $f(x) = a_0 + a_1x + a_2x^2 + ... a_nx^n$ 并且这里 $n - 1$ 是偶数。
我们把他按照奇偶性分开：
令 $g(x) = a_0 + a_2x^2 + a_4x^4 + ... a_{n - 2}x^{n - 2}$
令 $h(x) = a_1x + a_3x^3 + a_5x^5 + ... a_{n - 1}x^{n - 1}$
显然 $f(x) = g(x) + h(x)$ 这里 + 是按位相加。
再变：
$g(x)\rightarrow a_0 + a_2x + a_4x^2 + a_{n - 2}x^{\frac{n}{2} - 1}$
$h(x)\rightarrow a_1 + a_3x + a_5x^2 + a_{n - 1}x^{\frac{n}{2} - 1}$
所以现在 $f(x) = g(x ^ 2) + xh(x ^ 2)$，我们就可以带值找找规律了。

设当前带入 $\omega_n^k,k < \frac{n}{2}$ 分两种情况带入：

带入 $\omega_{n}^k$ 此时 $f(\omega_{n}^k) = g((\omega_{n}^k)^2) + \omega_{n}^kh((\omega_{n}^k)^2)$

$f(\omega_{n}^k) = g((\omega_{n}^k)^2) + \omega_{n}^kh((\omega_{n}^k)^2)\rightarrow f(\omega_{n}^k) = g(\omega_{n}^{2k}) + \omega_{n}^kh(\omega_{n}^{2k})\rightarrow f(\omega_{n}^k) = g(\omega_{\frac{n}{2}}^k) + \omega_{n}^kh(\omega_{\frac{n}{2}}^k)$。
由此得出来了，对于 $k < \frac{n}{2}$，$f(\omega_n^k)$ 的点值就是 $g(\omega_{\frac{n}{2}}^k) + \omega_n^kh(\omega_{\frac{n}{2}} ^ k)$ 也就是说，而 $g,h$ 都是子问题，且范围缩小了一倍，所以我们每次把原问题的范围分成两半来分别解决然后合并，这是典型的分治，总的层数是 $O(\log n)$ 级别的。

带入 $\omega_{n}^{k + \frac{n}{2}}$ 此时 $f(\omega_{n}^{k + \frac{n}{2}}) = g((\omega_{n}^{k + \frac{n}{2}}) ^ 2) + \omega_n^{k + \frac{n}{2}}h((\omega_n^{k + \frac{n}{2}}) ^ 2)$

因为 $(\omega_{n}^j)^k = \omega_n^{jk}$，原式得：$f(\omega_{n}^{k + \frac{n}{2}}) = g(\omega_{n}^{2k + n}) + \omega_n^{k + \frac{n}{2}}h(\omega_n^{2k + n})$
因为 $\omega_n^j = \omega_n^{j\bmod n}$，原式得：$f(\omega_{n}^{k + \frac{n}{2}}) = g(\omega_{n}^{2k}) + \omega_n^{k + \frac{n}{2}}h(\omega_n^{2k})$
因为 $\omega_{2n}^{2j} = \omega_n^j$，原式得：$f(\omega_n^{k + \frac{n}{2}}) = g(\omega_{\frac{n}{2}} ^ k) + \omega_n^{k + \frac{n}{2}} h(\omega_{\frac{n}{2}}^k)$
这里为了方便，可以根据 $\omega_n^{k + \frac{n}{2}} = -\omega_n^k$，原式得：$f(\omega_n^{k + \frac{n}{2}}) = g(\omega_{\frac{n}{2}} ^ k) - \omega_n^{k} h(\omega_{\frac{n}{2}}^k)$

因此我们就得到了关于计算 $f(\omega_n^i)$ 的式子，也就是每次把原序列折半分成奇数偶数两半，然后分别计算子问题的 $f(\omega_n^i)$ 最后合并，不难发现我们转换后的式子恰好满足 $n\rightarrow \frac{n}{2}$ 且 $k < \frac{n}{2}$ 所以成功划分成了子问题！

这里有些细节：

如果递归到某一层， $n = 1$ 那么直接 return，这个是因为递归到 $n = 1$ 的时候是计算 $f(w_1^0)$ 由于 $w_1 ^ 0$ 其实就是 $1 + 0i = 1$ 所以 $f(w_1 ^ 0) = $ 分治到 $i$ 时 $i$ 的系数，所以可以直接分治的时候记录系数，然后把系数当做 $f(w_1 ^ 0)$ 来计算。

时间复杂度不难发现是 $T(n) = 2T(\frac{n}{2}) + O(n) = O(n\log n)$

离散逆傅里叶变换 IDFT

假设我们已经进行了一次正变换 DFT，得到的 $N$ 个点值是序列 $\{y\},\forall i\in[1,|y|]$ 对于 $y_i = f(\omega_n^i)$ 设原来的函数是 $f(x) = a0 + a1x + a2x ^ 2 + a3x ^ 3 + ... a_{n - 1} x ^ {n - 1}$ 我们现在来构造一个新的多项式 $A$ 满足：

$A(x) = \sum_{0\leq i< n} y_ix^i$

我们再设 $b$ 序列，其中 $b_i = \omega_n^{-i}$ 那么我们来看看把 $\{b\}$ 序列带入 $A(x)$ 看看会产生什么样的点值序列。

$A(b_k) = \sum_{i = 0}^{n - 1} y_i{b_k}^i$

因为 $y_i = f(\omega_n ^ i) = \sum_{j = 0}^{n - 1} a_j({\omega_n ^ i}) ^ j$

所以原式得：$A(b_k) = \sum_{i = 0}^{n - 1} {b_k} ^ i \sum_{j = 0}^{ n- 1} a_j(\omega_n^{i})^j$

因为： $(\omega_n^i)^j = \omega_n ^ {ij}$

所以原式得：$A(b_k) = \sum_{i = 0}^{n - 1} {b_k} ^ i \sum_{j = 0}^{ n- 1} a_j\omega_n^{ij}$

因为 $b_k = \omega_n^{-k}$

所以原式得：$A(b_k) = \sum_{i = 0}^{n - 1} {\omega_n^{-k}} ^ i \sum_{j = 0}^{ n- 1} a_j\omega_n^{ij}\rightarrow \sum_{i = 0}^{n - 1} {\omega_n^{-ik}} \sum_{j = 0}^{ n- 1} a_j\omega_n^{ij}$

交换求和：
$A(b_k) = \sum_{i = 0}^{n - 1} \sum_{j = 0}^{ n- 1} {\omega_n^{-ik}}a_j\omega_n^{ij}$

因为：$\omega_n^i\omega_n^j = \omega_n^{i + j}$

所以原式得：$A(b_k) = \sum_{i = 0}^{n - 1} \sum_{j = 0}^{ n- 1} a_j\omega_n^{ij}\omega_n^{-ik}\rightarrow \sum_{i = 0}^{n - 1} \sum_{j = 0}^{ n- 1} a_j\omega_n^{i(j - k)}$

因为 $(\omega_n^i)^j = \omega_n^{ij}$

所以原式得：$A(b_k) = \sum_{i = 0}^{n - 1} \sum_{j = 0}^{ n- 1} a_j(\omega_n^{j -k}) ^ i$

交换求和得：
$A(b_k) = \sum_{j = 0}^{n - 1} a_j\sum_{i = 0}^{n - 1} (\omega_n^{j - k})^i$

此时记 $S(\omega_n^a) = \sum_{i = 0}^{n - 1} (\omega_n^{a}) ^ i$

如果 $a \bmod n = 0$ 那么显然 $\omega_n^a = 1$ 所以 $\sum_{i = 0}^{n - 1}(\omega_n^a) ^ i = n \times 1 = n$
否则考虑相邻两项作差：

等式两边同乘 $\omega_n^a$
原式得:$\omega_n ^ a S(\omega_n^a)=\sum_{i = 1}^{n}(\omega_n^{a}) ^ i$
两者作差得：$\omega_n^aS(\omega_n^a) - S(\omega_n ^ a) = \sum_{i = 1}^n (\omega_n^a)^i - \sum_{i = 0}^{n - 1}(\omega_n^a)^i$
左边提取公因式得：$(\omega_n^a - 1)S(\omega_n^a)= \sum_{i = 1}^n (\omega_n^a)^i - \sum_{i = 0}^{n - 1}(\omega_n^a)^i$
由于 $a\bmod n\ne 0$ 所以 $\omega_n^a\ne 1$ 所以 $(\omega_n^a - 1)\ne 0$ 所以可以等式两边同时除以 $(\omega_n^a - 1)$
因此，原式得：$S(\omega_n^a) = \frac{\sum_{i = 1}^n (\omega_n^a)^i - \sum_{i = 0}^{n - 1}(\omega_n^a)^i}{\omega_n^a - 1}$
不难发现 $\sum_{i = 1}^n (\omega_n^a)^i - \sum_{i = 0}^{n - 1}(\omega_n^a)^i = (\omega_n^a)^n - (\omega_n^a)^0$
前者可以写成 $(\omega_n^n)^a = (\omega_n^0)^n = 1$ 后者是 $0$ 次方也是 $1$
所以可以得到对于 $a\bmod n \ne 0$ 的所有 $a$ 满足 $S(\omega_n^a) = 0$

所以最后可以归纳成：对于 $a < n$ 的所有 $\omega_n^a$ 满足：也就是说

$ S\left(\omega_n^a\right)=
\begin{cases}
n,&a=0\
0,&a\neq 0
\end{cases}
$

好的，现在我们把这个公式带到 $A(b_k) = \sum_{j = 0}^{n - 1} a_j\sum_{i = 0}^{n - 1} (\omega_n^{j - k})^i$

原式得：$\sum_{j = 0}^{n - 1} a_j S(\omega_n^{j - k})$

显然这个式子，只有在 $j = k$ 的时候有值，且值为 $n$ 所以这个式子的值就是 $na_k$

也就是说，给定点 $b_i = \omega_n ^ {-i}$ 它的点值序列表示就是 $\{(b_0,a_0n),(b_1,a_1n),....,(b_{n - 1},a_{n - 1}n)\}$

所以我们就得到了如何在求出 $f,g$ 的点值序列后如何还原出来新的多项式的系数，当然，因为最后还原出来的点值是 $a_in$ 所以要除以 $n$。

不难发现，我们这个过程和 DFT 几乎一模一样，无非就是将 $n$ 个单位根 $\omega_n^i\rightarrow \omega_n^{-i}$ 这两个式子是关于实轴（x轴）对称的，所以只需要把纵坐标乘上一个 $-1$ 即可。

但是我们发现直接写分治，我们空间复杂度应该也是 $O(n \log n)$ 的，因为需要对每个点开一个 $O(len)$ 级别的数组，所以要引进一种空间上的优化：蝶形优化。

蝶形优化

对于当前分治区间 $[l_i,r_i]$ 设当前长度是 $n_i = r_i - l_i + 1$

我们要计算 $f(\omega_{n_i}^{0,1,...n_i - 1})$ 需要 $g(\omega_{\frac{n_i}{2}} ^ {0,1,2,...\frac{n_i}{2} - 1})$ 和 $h(\omega_{\frac{n_i}{2}} ^ {0,1,2,...\frac{n_i}{2} - 1})$

我们可以在递归开始前，就把 $g,h$ 给放到 $f$ 的第 $[0 , \frac{n_i}{2} - 1]$ 的位置和 $[\frac{n_i}{2},n_i - 1]$ 的位置上，这个位置是个指针，此时修改 $g$ 和 $h$ 也会在对应的 $f$ 上进行修改，（类似长剖一样）所以在我们更新 $f$ 的时候不能直接拿 $g$ 和 $h$ 来改 $f$ 而是新开一个辅助数组 $G$ 最后再把 $G$ 赋值给 $f$ 不然直接改 $f$ 的值，对应的 $g$ 和 $h$ 的值也会进行修改，就不对了。这样空间就被优化成 $O(n)$ 级别了。

最后，一些实现上的小细节：

我们发现当 $n = 1$ 的时候 $\omega_n^{n - 1} = 1$ 的，所以可以直接把一开始的点值放进去，因为此时点值的 $f$ 值就是自己不会变。所以一开始的时候，就直接输入 $A,B$ 的实部就可以了。
一定要注意，$f(i) = a + bi$ 中 $a + bi$ 中的 $(a,b)$ 是个整体，代表一个复数，而不是当 $x = a$ 的时候对应的值是 $b$。
分治的时候，一定不要直接修改 $f$ 的值，因为使用蝶形优化后，$g$ 和 $h$ 和 $f$ 公用一个数组。
因为我们刚才说的，每次分治的时候都需要能均分奇偶数，所以要求所有大于 $1$ 的分治区间长度是 $2$ 的倍数，自然想到构造一个 $2 ^ k$ 的数列，如果一开始不够就补 $0$ 即可，然后找到最大的 $k$ 满足 $2^k \ge n + m + 1$
由于我们是进行的实数运算，会有精度问题，所以最后要输出四舍五入后的值。

代码：

#include <iostream>
#include <cmath>
#include <algorithm>constexpr double Pi = acos(-1);
constexpr int N = 4e6 + 5;
#define int long long
struct complex{double x,y;inline complex operator + (const complex & a){return {a.x + x,a.y + y};}inline complex operator - (const complex & a){return {x - a.x,y - a.y};}inline complex operator * (const complex & a){// (x + yi) * (a.x + a.yi)// a.x * x + xa.yi + a.xyi + a.yi*yi// a.x * x + xa.yi + a.xyi - a.y*yreturn {a.x * x - a.y * y,x * a.y + a.x * y};}
}A[N],B[N],g[N];
int n , m;
inline void fft(complex * f,int len,int rev){if(len == 1) return;for(int i = 0; i < len; ++ i){g[i] = f[i];}complex * ls = f,* rs = f + len / 2;for(int i = 0; i < len / 2; ++ i){ls[i] = g[i << 1];rs[i] = g[i << 1 | 1];}fft(ls,len / 2,rev);fft(rs,len / 2,rev);complex w = {cos(2 * Pi / len),sin(2 * Pi / len) * 1.0 * rev},now;now.x = 1;now.y = 0;for(int i = 0; i < len / 2; ++ i){g[i] = ls[i] + now * rs[i];g[i + len / 2] = ls[i] - now * rs[i];now = now * w;}for(int i = 0; i < len; ++ i){f[i] = g[i];}
}
signed main(){std::ios::sync_with_stdio(false);std::cin.tie(nullptr);std::cin >> n >> m;for(int i = 0; i <= n; ++ i){std::cin >> A[i].x;}for(int i = 0; i <= m; ++ i){std::cin >> B[i].x;}int len = 1;for(int k = n + m; len <= k; len <<= 1);fft(A,len,1);fft(B,len,1);for(int i = 0; i <= len; ++ i){A[i] = A[i] * B[i];}fft(A,len,-1);for(int i = 0; i <= n + m; ++ i){std::cout << (int)(A[i].x / len + 0.5) << ' ';}return 0;
}

下面介绍迭代法实现快速傅立叶变换

我们发现，其实 $0 $ 到 $n - 1$ 的最后位置是已经确定好的了，比如 $n = 8$ 序列 $\{0,1,2,3,4,5,6,7\}$ 分治的过程如下：

原序列： $\{0,1,2,3,4,5,6,7\}$

一次变换： $\{0,2,4,6\},\{1,3,5,7\}$

二次变换： $\{0,4\},\{2,6\},\{1,5\},\{3,7\}$

三次变换： $\{0\},\{4\},\{2\},\{6\},\{1\},\{5\},\{3\},\{7\}$

最后序列： $\{0,4,2,6,1,5,3,7\}$

我们发现了什么规律，没错，如果设 $2^M = n + 1$ 那么设如果二进制一共 $M$ 位，也就是 $[0,n - 1]$ 都表示成 $M$ 位多项式，如果位数不够就前面补 $0$，那么最终的位置就是 $i$ 二进制下取反再转成十进制表示的位置，我们把这样的操作叫位逆序置换。

如何求解位逆序置换：
显然有一个 $O(n\log n)$ 的做法，对每个数暴力分解二进制然后翻转，单个数时间复杂度是 $O(\log n)$ 的，总的就是 $O(n\log n)$。
但是其实可以 $O(n)$ 从小到大地求出位逆序置换。
我们称位逆序置换后的数组为 $R$
首先 $R_0 = 0$
假设我们已经知道了 $R(\lfloor \frac{n}{2}\rfloor)$ 我们能否求出 $R(n)$ 呢？我们发现，其实 $R(n)$ 就是 $R(\lfloor \frac{n}{2}\rfloor)$ 二进制下右移一下，但是这充分了吗？
正确规律应该是，如果一个数右移了一位，且如果 $x$ 末位是 $1$ 那么就把最高位变成 $1$，这里因为是右移，所以之前最高位一定是 $0$，所以可以直接加上 $2^{M - 1}$
所以我们可以这样计算一个 $R(n)$，$R(n) = \lfloor\frac{R(\lfloor\frac{n}{2}\rfloor)}{2}\rfloor$ 也就是 C++ 中的右移一位，然后判断一下此时的 $R(n)$ 末位是否是 $1$，如果是 $1$ 就让 $R(n)\rightarrow R(n) + 2 ^ {M - 1}$
使用 C++ 中的位运算可以让单次操作复杂度降至 $O(1)$ 总的时间复杂度为 $O(n)$
值得注意的是，最后我们求出来了 $R$ 数组，并且要 $i$ 和 $R$ 交换位置，实际上我们只需要在 $i < R_i$ 的时候交换就可以了，不然在 $R_i$ 再进行交换就相互抵消了（我的理解）。

所以我们可以直接把每个位置的最后位置求出来，然后迭代法合并结果即可，具体就是先从小到大枚举当前分治区间的长度（这里不是真的分治，是类比之前的分治做法），第二轮枚举在这个区间长度的前提下所有满足分治区间长度是当前枚举区间长度的区间左端点，最后一维则是枚举每个当前长度左端点的对应区间，复杂度也是 $O(n\log n)$ 的，但是会比递归要快一些。

迭代版本代码，其中 reverse 函数就是实现求解 $R$ 的过程。

#include <iostream>
#include <cmath>using ll = long long;constexpr int N = 4e6 + 5;
constexpr double Pi = acos(-1);class complex{public:double x,y;complex operator + (const complex & a){return (complex){x + a.x,y + a.y};}complex operator - (const complex & a){return (complex){x - a.x,y - a.y};}complex operator * (const complex & a){// (x + yi) * (a.x + a.yi)// a.x * x + x * a.yi + a.x * yi + a.yi * yi// i ^ 2 = - 1// a.x * x  - a.y * y + (a.y * x + a.x * y) ireturn (complex){a.x * x - a.y * y,a.y * x + a.x * y};}
}a[N],b[N];
int n,m,len = 1,R[N];
inline void reverse(complex f[],int len){R[0] = 0;for(int i = 1; i < len; ++ i){R[i] = R[i >> 1] >> 1;if(i & 1){R[i] |= len >> 1;}}for(int i = 0; i < len; ++ i){if(i < R[i]){std::swap(f[i],f[R[i]]);}}
}
inline void fft(complex f[],int len,int rev){reverse(f,len);for(int L = 2; L <= len; /*下标是 0 ~ n - 1 长度是 n*/ L <<= 1){complex w_1 = {cos(2 * Pi / L),sin(2 * Pi / L) * rev};for(int l = 0; l < len; l += L){complex now = {1,0};for(int i = l; i < l + L / 2; ++ i,now = now * w_1){complex g = f[i],h = now * f[i + L / 2];f[i] = g + h;f[i + L / 2] = g - h;}}}
}
signed main(){std::ios::sync_with_stdio(false);std::cin.tie(nullptr);std::cin >> n >> m;for(int i = 0; i <= n; ++ i){std::cin >> a[i].x ;// 读入实部 }for(int i = 0; i <= m; ++ i){std::cin >> b[i].x;}for(int k = n + m; len <= k; len <<= 1);fft(a,len,1);fft(b,len,1);for(int i = 0; i <= len; ++ i){a[i] = a[i] * b[i];}fft(a,len,-1);for(int i = 0; i <= n + m; ++ i){std::cout << (ll)(a[i].x * 1.0 / len + 0.5) << ' ';}return 0;
}

接下来讲解另一种多项式乘法的方法：

快速数论变换（NTT）

我们先来回顾一下 FFT，FFT 利用了比较特殊，且有周期性的单位根 $\omega_n^i$ 从而快速求出 $n$ 个点值，但是 FFT 运用的是复数，运算多是浮点数运算，具有精度问题，那么有没有一种算法，能有效规避这种精度而且具有和单位根一样良好的性质呢？

答曰：有的兄弟，有的，没错这就是 NTT

首先我们需要知道原根：

设一个质数 $p$ 和一个数 $a$，满足 $\gcd(a,p) = 1$，且存在 $x$ 满足 $a^x \bmod p = 1$ 那么我们就称 $x$ 是，称作 $a$ 模 $m$ 的阶，记作 $\delta_m(a)$ 或 $\operatorname{ord}_m(a)$

原根：

设 $m \in \mathbf{N}^{*}， g\in \mathbf{Z}$ 若
$(g,m)=1$，且
$\delta_m(g)=\varphi(m)$，则称
$g$ 为模
$m$ 的原根。

即 $g$ 满足
$\delta_m(g) = \left| \mathbf{Z}_m^* \right| = \varphi(m)$

最重要的性质之一：

若一个数 $m$ 有原根 $g$，则 $g,g^2,\ldots,g^{\varphi(m)}$ 构成模 $m$ 的简化剩余系。

特别地，当 $m$ 是质数时，有 $g^i\bmod m，0 < i < m$ 的结果两两不同，这也是 NTT 运用原根的核心部分。

现在我们来考虑把原根和单位根联系起来：

周期性：
- 单位根 $\omega$ 满足 $\omega_n^i = \omega_n^{i + n}$
- 原根 $g$ 满足 $g^i = g^{i + p - 1},p$ 是质数，这是因为费马小定理 $g^{i + p - 1} = g^i g^{p - 1}$，而 $g^{p - 1} \mod p = 1$
- 所以不难想到我们应该让每 $p - 1$ 作为一个周期，类比单位根，我们应该以 $\frac{p - 1}{n}$ 作为 $\omega_n^1$
合法性（乱取的）：
- 由于我们确定 $n$ 次多项式需要 $n + 1$ 个点值，且满足横坐标两两不同，而原根则满足这一点，具体的，因为 $1 < \frac{p - 1}{n}$ 且 $n\times (\frac{p - 1}{n}) = p - 1 < p$ 所以自然满足这 $n$ 个单位根两两不同。
可逆性（乱取的）：
- 注意到 FFT 的逆运算只需要把单位根变成关于实轴对称点，也就是 $\omega_n^i\rightarrow \omega_n^{-i}$ 其实说白了就是给做了个 $-1$ 次方，因为 $p$ 是质数，所以 NTT 中的 $\omega_n^i$ 也是有模意义下的乘法逆元的，最后除以 $n$ 也只需要乘上 $n$ 的逆元即可。

那么此时 NTT 中的“单位根”就有和 FFT 中的复数单位根有着相同的作用了。

由于只需要实数（整数）的运算，时间效率明显快了一些。

NTT 迭代版代码：

#include <iostream>using ll = long long;constexpr int N = 4e6 + 5;
constexpr ll mod = 998244353;ll a[N],b[N];
int R[N],n,m,len = 1;
inline void reverse(ll f[],int len){R[0] = 0;for(int i = 1; i < len; ++ i){R[i] = R[i >> 1] >> 1;if(i & 1) R[i] |= len >> 1;}for(int i = 0; i < len; ++ i){if(i < R[i]){std::swap(f[i],f[R[i]]);}}
}
inline ll ksm(ll x,int y){ll ans = 1;for(;y;y >>= 1,(x *= x) %= mod){if(y & 1) (ans *= x) %= mod;}return ans;
}
inline void NTT(ll f[],int len,int rev){reverse(f,len);for(int L = 2; L <= len; L <<= 1){ll w_1 = ksm(3,(mod - 1) / L);if(rev == -1) w_1 = ksm(w_1,mod - 2);for(int l = 0; l < len; l += L){ll now = 1;for(int i = l; i < l + L / 2; ++ i,(now *= w_1) %= mod){ll g = f[i],h = f[i + L / 2] * now % mod;f[i] = (g + h) % mod;f[i + L / 2] = (g - h + mod) % mod;}}}
}
int main(){std::ios::sync_with_stdio(false);std::cin.tie(nullptr);std::cin >> n >> m;for(int i = 0; i <= n; ++ i){std::cin >> a[i];}for(int i = 0; i <= m; ++ i){std::cin >> b[i];}for(int k = n + m; len <= k; len <<= 1);NTT(a,len,1);NTT(b,len,1);for(int i = 0; i < len; ++ i){(a[i] *= b[i]) %= mod;}NTT(a,len,-1);ll inv = ksm(len,mod - 2);for(int i = 0; i <= n + m; ++ i){(a[i] *= inv) %= mod;std::cout << a[i] << ' ';}return 0;
}

分治 FFT（CDQ + NTT）

前置只是：

CDQ 分治思想，快速数论变换（NTT）

题意描述：

给定序列 $g_{1,...,n - 1}$ 求序列 $f_{0,...,n - 1}$，其中 $f_i = \sum_{j = 1} ^ i g_j f_{i - j}$

首先，我们发现 $f$ 的计算方式比较特殊，想计算 $f_i$ 的时候，必须要知道 $f_{0,...i - 1}$，那么不难想到我们可以类似 CDQ 分治那样去做，就是先计算分治区间左边一半的 $f$，然后去更新右边的 $f$，再递归右边，而 CDQ 分治正确的原因就是左边的区间已经计算好了答案，不会存在左侧没计算完就更新右边的情况。

设当前我们分治区间是 $[l,r]$，设中点是 $mid$，我们现在应该计算 $f_{l,...mid}$ 和 $g$ 的前 $mid - l + 1$ 项对 $[mid + 1,r]$ 的 $f$ 的影响，我们发现，此时已经转换成了对 $f_{l,...mid}$ 和 $g$ 的前 $mid - l + 1$ 项做卷积，这里要认为 $f_l$ 是第一个多项式的第一项的系数，然后这样卷积后的第 $i$ 项，实际上是 $l + i$ 的位置，因为刚开始我们让 $f$ 往前移了 $l$ 位，乘法的时候使用 NTT 即可。

根据主定理，时间复杂度是 $T(n) = 2T(n / 2) + O(n\log n) = O(n\log ^ 2n)$，可以通过。

代码：

#include <iostream>using ll = long long;constexpr int N = 4e5 + 5;
constexpr ll mod = 998244353;int n,len = 0,R[N];
ll g[N],f[N],a[N],b[N];
inline ll ksm(ll x,ll y){ll ans = 1;for(;y;y >>= 1,(x *= x) %= mod){if(y & 1){(ans *= x) %= mod;}}return ans;
}
inline void reverse(ll *f,int len){R[0] = 0;for(int i = 1; i < len; ++ i){R[i] = R[i >> 1] >> 1;if(i & 1){R[i] |= len >> 1;}}for(int i = 0; i < len; ++ i){if(i < R[i]) std::swap(f[i],f[R[i]]);}
}
inline void NTT(ll *f,int len,int rev){reverse(f,len);for(int L = 2; L <= len; L <<= 1){ll w_1 = ksm(3,(mod - 1) / L);if(rev == -1){w_1 = ksm(w_1,mod - 2);}for(int l = 0; l < len; l += L){ll w_now = 1;for(int i = l; i < l + L / 2; ++ i){ll g = f[i],h = w_now * f[i + L / 2] % mod;f[i] = (g + h) % mod;f[i + L / 2] = (g - h + mod) % mod;(w_now *= w_1) %= mod;}}}if(rev == -1){for(int i = 0; i < len; ++ i){f[i] = f[i] * ksm(len,mod - 2) % mod;}}
}
inline void solve(int l,int r){if(l == r){if(!l) f[l] = 1;return;}int mid = (l + r) >> 1;solve(l,mid);for(int i = l; i <= mid; ++ i){a[i - l] = f[i];}for(int i = 1; i <= r - l + 1; ++ i){b[i] = g[i];}len = 1;for(;len <= r - l + 1; len <<= 1);for(int i = mid - l + 1; i < len; ++ i){a[i] = 0;}for(int i = r - l + 2; i < len; ++ i){b[i] = 0;}NTT(a,len,1);NTT(b,len,1);for(int i = 0; i < len; ++ i){(a[i] *= b[i]) %= mod;}NTT(a,len,-1);for(int i = mid + 1; i <= r; ++ i){(f[i] += a[i - l]) %= mod;}solve(mid + 1,r);
}
int main(){std::ios::sync_with_stdio(false);std::cin.tie(nullptr);std::cin >> n;for(int i = 1; i < n; ++ i){std::cin >> g[i];}solve(0,n - 1);for(int i = 0; i < n; ++ i){std::cout << f[i] << ' ';}return 0;
}

多项式求逆

前置知识：

快速数论变换（NTT），模意义下的乘法逆元。

题意描述：

给定 $n - 1$ 次多项式 $f(x) = \sum_{i = 0}^{n - 1} a_ix^i$，求一个 $n - 1$ 次多项式 $g$，满足 $f\times g \equiv 1 \pmod {x ^ n}$

首先，一个 $n$ 次多项式 $f$，对 $x ^ k$ 次方取模的结果如下：

若 $k > n$ 那么取模后 $f$ 不变。
否则 $f\rightarrow \sum_{i = 0}^{k - 1} a_i x ^ k$ 也就是保留 $f$ 的前 $k$ 项。

那么由于本题是 $n - 1$ 次多项式，而要求对 $x ^ n$ 取模，所以等价于找到一个多项式 $g$，满足 $g$ 只有常数项是 $1$，其他项系数都是 $0$

考虑分治处理这个问题

首先不妨设 $n = 2 ^ k,k\in N$

下面分两种情况：

当 $n = 1$ 的时候，就是求解 $f_0\times g_0\equiv 1\pmod {x ^ 1}$，显然此时 $g_0 = {f_0}^{-1}$，求 $f_0$ 的逆元作为 $g_0$ 即可。
当 $n > 1$ 的时候：
- 假设我们已经求出来了 $f$ 在模 $x^{\frac{n}{2}}$ 的时候符合要求的 $g$，记为 $g'$。
- 那么自然有 $g'f \equiv 1\pmod{x ^ {\frac{n}{2}}}$
- 即：$g'f - 1\equiv 0\pmod{x ^ {\frac{n}{2}}}$
- 由于 $x\equiv 0\pmod{p}$ 那么有 $x^2 \equiv 0\pmod {p ^ 2}$
  - 证明：因为 $x\equiv 0\pmod{p}$ 所以 $x\bmod p = 0$ 即，$x = kp,k\in Z$，那么 $x ^ 2 = k ^ 2p ^ 2 = (kp) ^ 2$，因为 $k\in Z$ 所以 $k ^ 2\in Z$，所以 $x ^ 2\bmod p = 0$，所以 $x ^ 2\equiv 0\pmod{p ^ 2}$。
- 那么我们把左右两边同时平方，根据完全平方公式得：$(g'f) ^ 2 + 1 - 2g'f\equiv 0\pmod{x ^ n}$，即：$(g'f)^2 - 2g'f \equiv -1\pmod{x ^ n}$
- 两边乘 $-1$ 得：$2g'f - (g'f) ^ 2\equiv 1\pmod{x ^ n}$
- 设 $gf\equiv 1\pmod{x ^ n}$
- 那么联立上方两个式子，得到：$(g'f) ^ 2 - 2g'f\equiv gf\pmod{x ^ n}$
- 两边同时除以 $f$ 得：$f{g'} ^ 2 - 2g'\equiv g\pmod{x ^ n}$，那么根据这个公式就可以求解 $g$ 了。

由于我们计算模 $x ^ n$ 的符合的 $g$ 时只需要知道模 $x ^ {\frac{n}{2}}$ 的符合的 $g'$

所以我们时间复杂度根据主定理是 $T(n) = T(n / 2) + O(n\log n) = O(n\log n)$

一个细节就是，递归的时候长度应该要设成 $2 \times n$ 而不是 $n$，因为 NTT 的时候需要二的幂次，而 $f$ 是 $n$ 项的，$g'$ 是 $\frac{n}{2}$ 项的，取二的整次幂就是 $2\times n$

思考：如果求一个 $n - 1$ 次多项式 $f(x)$ 此时不是对 $x ^ n$ 取模，而是对 $m < n$ 的 $x ^ m$ 取模呢？

我们发现此时只有 $f(x)$ 的前 $m$ 项式有用的，所以直接等价于求解一个 $m - 1$ 次多项式，对 $x ^ m$ 取模了。

代码：

#include <iostream>constexpr int N = 4e5 + 5;
using ll = long long;
constexpr ll mod = 998244353;
ll ksm(ll x,ll y){ll ans = 1;while(y){if(y & 1) (ans *= x) %= mod; (x *= x) %= mod;y >>= 1;}return ans;
}
ll inv(ll x){return ksm(x,mod - 2);
}
int R[N * 2];
inline void reverse(ll f[],int n){R[0] = 0;for(int i = 1; i < n; ++ i){R[i] = R[i >> 1] >> 1;if(i & 1){R[i] |= n >> 1;}}for(int i = 0; i < n; ++ i){if(i < R[i]){std::swap(f[i],f[R[i]]);}}
}
inline void NTT(ll f[],int n,int rev){reverse(f,n);for(int L = 2; L <= n; L <<= 1){ll w_1 = ksm(3,(mod - 1) / L);if(rev == -1){w_1 = inv(w_1);}for(int l = 0; l < n; l += L){ll w_now = 1;for(int i = l; i < l + L / 2; ++ i){ll g = f[i],h = w_now * f[i + L / 2] % mod;f[i] = (g + h) % mod;f[i + L / 2] = (g - h + mod) % mod;(w_now *= w_1) %= mod; }}}if(rev == -1){ll inv_n = inv(n);for(int i = 0; i < n; ++ i){(f[i] *= inv_n) %= mod;}}
}
ll c[N],d[N];
inline void solve(ll a[],ll b[],int n){if(n == 1){b[0] = inv(a[0]);return;}for(int i = 0; i <= n * 2; ++ i){c[i] = d[i] = 0;}solve(a,c,n >> 1);for(int i = 0; i < n; ++ i){d[i] = a[i];}NTT(d,n << 1,1);NTT(c,n << 1,1);for(int i = 0; i < n * 2; ++ i){b[i] = (2 * c[i] - d[i] * c[i] % mod * c[i] % mod + mod) % mod;}NTT(b,n << 1,-1);for(int i = 0; i < n; ++ i){b[i + n] = 0;}
}
int n,len = 0;
ll a[N],b[N];
int main(){std::ios::sync_with_stdio(false);std::cin.tie(nullptr);std::cin >> n;for(int i = 0; i < n; ++ i){std::cin >> a[i];}for(len = 1; len <= n; len <<= 1);solve(a,b,len);for(int i = 0; i < n; ++ i){std::cout << b[i] << ' ';}return 0;
}

多项式除法

定义：

$f = gq + r$

其中，如果 $f$ 是一个 $n$ 次多项式，$g$ 是 $m$ 次多项式，那么 $q$ 是 $n - m$ 次多项式，$h$ 是小于 $m$ 次多项式

我们令 $f' = rev(f)$ 即把 $f$ 的系数翻转，$f'(x) = \sum_{i = 0}^n a_ix^{n - i}$，其实 $f'(x) = f(\frac{1}{x})x ^ n$

那么有 $f(x) = g(x)q(x) + r(x)$

把 $\frac{1}{x}$ 带入得到：$f(\frac{1}{x}) = g(\frac{1}{x})q(\frac{1}{x}) + r(\frac{1}{x})$

等式两边同时乘上 $x ^ n$ 得到：$f(\frac{1}{x})x^n = g(\frac{1}{x})x ^ mq(\frac{1}{x})x ^ {n - m} + r(\frac{1}{x})x ^ n$

$f' = g'q' + r'x^{n - m + 1}$，这里因为 $deg_r < m$ 所以最大是 $m - 1$，那么剩余的幂次就是 $n - (m - 1) = n - m + 1$

那么移项得：$f' - g'q' = r'x^{n - m + 1}$

注意到此时如果我们对等式两边同时对一个数 $x$ 取模，等式两边是同余的，那么我们让等式两边同时对 $x ^ {n - m + 1}$ 取模，由于等式右边最低次幂是 $n - m + 1$ 所以等式右边对 $x ^ {n - m + 1}$ 取模后，右边应该是 $0$，所以 $f' - g'q'\equiv 0\pmod{x ^ {n - m + 1}}$

移项得：$f'\equiv g'q'\pmod{x ^ {n - m + 1}}$

同时除以 $g'$ 得：

$\frac{f'}{g'}\equiv q'\pmod{x ^ {n - m + 1}}$

注意到 $q'$ 的次数是 $n - m$ 所以可以直接把 $q'$ 的系数再翻转回去求出 $q$

那么知道 $q$，很容易也能求出 $r$ 了，也就是 $r = f - gq$

时间复杂度是求逆元的复杂度为 $O(n\log n)$

代码：

#include <iostream>
#include <algorithm> #define int long long
constexpr int N = 6e5 + 5,mod = 998244353;
int n,m,f[N],g[N],rev_f[N],rev_g[N],inv_g[N],q[N],r[N];
inline int ksm(int x,int y){int ans = 1;for(;y; y >>= 1,(x *= x) %= mod){if(y & 1) (ans *= x) %= mod;}return ans;
}
struct poly{int R[N],c[N * 2],d[N * 2];inline void reverse(int f[],int n){R[0] = 0;for(int i = 1; i < n; ++ i){R[i] = R[i >> 1] >> 1;if(i & 1){R[i] |= n >> 1;}}for(int i = 0; i < n; ++ i){if(i < R[i]){std::swap(f[i],f[R[i]]);}}}inline void NTT(int f[],int n,int rev){reverse(f,n);for(int L = 2; L <= n; L <<= 1){int w_1 = ksm(3,(mod - 1) / L);if(rev == -1) w_1 = ksm(w_1,mod - 2); for(int l = 0; l < n; l += L){int w_now = 1;for(int i = l; i < l + L / 2; ++ i,(w_now *= w_1) %= mod){int g = f[i],h = w_now * f[i + L / 2] % mod;f[i] = (g + h) % mod;f[i + L / 2] = (g - h + mod) % mod;}}}if(rev == -1){int inv_n = ksm(n,mod - 2);for(int i = 0; i < n; ++ i){(f[i] *= inv_n) %= mod;}}}inline void inv(int a[],int b[],int len){if(len == 1){b[0] = ksm(a[0],mod - 2);return;}for(int i = 0; i <= (len << 1); ++ i){c[i] = d[i] = 0;}inv(a,c,(len + 1) >> 1);for(int i = 0; i < len; ++ i){d[i] = a[i];}NTT(c,len << 1,1);NTT(d,len << 1,1);for(int i = 0; i < (len << 1); ++ i){b[i] = (c[i] * 2 % mod - c[i] * c[i] % mod * d[i] % mod + mod) % mod;}NTT(b,len << 1,-1);for(int i = 0; i < len; ++ i){b[i + len] = 0;}}inline void Inv(int f[],int g[],int n){int len = 1;for(;len <= n; len <<= 1);inv(f,g,len);}inline void mul(int f[],int g[],int n){int len = 1;for(;len <= n; len <<= 1);NTT(f,len,1);NTT(g,len,1);for(int i = 0; i < len; ++ i){(f[i] *= g[i]) %= mod;}NTT(f,len,-1);}
}P;
signed main(){std::ios::sync_with_stdio(false);std::cin.tie(nullptr);std::cin >> n >> m;for(int i = 0; i <= n; ++ i){std::cin >> f[i];rev_f[n - i] = f[i];}for(int i = 0; i <= m; ++ i){std::cin >> g[i];rev_g[m - i] = g[i];}for(int i = n - m + 1; i <= m; ++ i){rev_g[i] = 0;}P.Inv(rev_g,inv_g,n - m + 1);P.mul(rev_f,inv_g,(n << 1) + m);for(int i = 0; i <= n - m; ++ i){q[i] = rev_f[n - m - i];std::cout << q[i] << ' ';}std::cout << '\n';P.mul(g,q,n);for(int i = 0; i < m; ++ i){r[i] = (f[i] - g[i] + mod) % mod;std::cout << r[i] << ' ';}return 0;
}

FWT 快速沃尔什变换

FWT 通常用于处理位运算的卷积运算，形如 $C_i = \sum_{i = j + k} A_j . B_k$ 其中，$.$ 分别代表 $\mid,\text{\&},\oplus$ 也就是或，与和异或。

FWT 具有线性性，即 $FWT(A + B) = FWT(A) + FWT(B)$，还有 $FWT(\lambda A) = \lambda FWT(A)$

下面来看这个题：

分别求：

$C_i = \sum_{j\mid k = i} A_j B_k\\ C_i = \sum_{j\text{\&} k = i} A_j B_k\\ C_i = \sum_{j\oplus k = i} A_j B_k$

我们分别来看看如何求解：

或（or）卷积：

类比 FFT 的方法，我们让 $A\rightarrow FWT_A,B\rightarrow FWT_B$，然后令 $FWT_C = FWT_A\times FWT_B$，然后让 $FWT_C\rightarrow C$ 即可得到 $C$

对于或卷积，令 ${FWT_A}_i = \sum_{j\subseteq i} A_j$，同理设出 $FWT_B$，下面证明 $FWT_C = FWT_A\times FWT_B$

证明：

${FWT_C}_k = {FWT_A}_k \times {FWT_B}_k$

${FWT_A}_k \times {FWT_B}_k = (\sum_{p\subseteq k} A_p)(\sum_{q\subseteq k } B _q)$
$p\subseteq k,q\subseteq k\Leftrightarrow p\mid q\subseteq k$
${FWT_A}_k \times {FWT_B}_k = (\sum_{p\subseteq k} A_p)(\sum_{q\subseteq k } B _q) = \sum_{p\mid q\subseteq k} A_p B_q = {FWT_C}_k$

根据这个式子我们可以分治解决，设当前分治区间前一半是 $A_0$，后一半是 $A_1$，那么 $FWT$ 当前区间等于 $\text{merge}(FWT_{A_0},FWT_{A_0} + FWT_{A_1})$ 这里 + 是按位相加，merge 类似字符串拼接，把前后拼起来组成 $FWT$ 序列。

说一下为啥，注意到对于 $A_0$ 的第 $i$ 个数，他对应的 $A_1$ 的第 $i$ 个数，在二进制表示下只是比 $A_0$ 最高位多了一个 $1$，所以 $A_0$ 的第 $i$ 个数也是 $A_1$ 的第 $i$ 个数的子集，也就是 $A_0\subseteq A_1$ 所以要把 $FWT_{A_0}$ 的贡献也要加到 $FWT_{A_1}$ 的对应位置上。

逆变换：

考虑如何根据 $FWT_A$ 来还原 $A$ 的值，注意到 $\text{merge}$ 的时候是 $\text{merge}(FWT_{A_0},FWT_{A_0} + FWT_{A_1})$，那么我们要还原就是把 $A_1$ 的位置减去之前 $FWT_{A_0}$ 的贡献。

于是可以和 FFT 一样去迭代计算即可，或卷积的代码：

inline void FWT_or(int *f,int rev){for(int L = 2; L <= (1 << n); L <<= 1){for(int l = 0; l < (1 << n); l += L){for(int i = l; i < l + L / 2; ++ i){(f[i + L / 2] += f[i] * rev % mod) %= mod;}}}
}

接下来讲与卷积：

类比或卷积，我们设 ${FWT_A}_i = \sum_{i\subseteq j} A_j$，也就是 $i$ 的超集的和。

下面证明 ${FWT_C}_k = {FWT_A}_k \times {FWT_B}_k$

证明：

${FWT_A}_k\times {FWT_B}_k = (\sum_{k\subseteq p} A_p)(\sum_{k\subseteq q} B_q)$
$k\subseteq p,k\subseteq q\Leftrightarrow k\subseteq p\text{\&} q$
所以 ${FWT_A}_k\times {FWT_B}_k = (\sum_{k\subseteq p} A_p)(\sum_{k\subseteq q} B_q) = \sum_{k\subseteq p\text{\&} q} A_p B_q = {FWT_C}_k$

那么类比刚刚的 $\text{merge}$，得到现在的 $\text{merge}$ 就是 $FWT(A) = \text{merge}(FWT_{A_0} + FWT_{A_1},FWT_{A_1})$

这个也很好解释，就是对于 $A_0$ 的第 $i$ 个数，他和 $A_1$ 的第 $i$ 个数的唯一差别就是最高位没有 $1$（这个同上面的或卷积），所以 $A_1$ 的第 $i$ 个数包含的超集 $A_0$ 的第 $i$ 个数也包含，所以有这个 $\text{merge}$ 的过程。

与卷积的代码：

inline void FWT_and(int *f,int rev){for(int L = 2; L <= (1 << n); L <<= 1){for(int l = 0; l < (1 << n); l += L){for(int i = l; i < l + L / 2; ++ i){(f[i] += f[i + L / 2] * rev % mod) %= mod;}}}
}

异或卷积：

这个卷积是最难的也是最妙的，首先定义 $x\circ y = popcount(x\text{\&} y)\bmod 2$，定义 ${FWT_A}_i = \sum_{j\circ i = 0} A_j - \sum_{k\circ i = 1} A_k$

$FWT[A]_i=\sum_{i\circ j=0}A_j-\sum_{i\circ j=1}A_j$。我们来证一下
$FWT[C] = FWT[A] \cdot FWT[B]$ 的正确性：

$\begin{aligned} FWT[A]_iFWT[B]_i&=\left(\sum_{i\circ j=0}A_j \sum_{i\circ j=1}A_j\right)\left(\sum_{i\circ k=0}B_k-\sum_{i\circ k=1}B_k\right) \\ &=\left(\sum_{i\circ j=0}A_j\sum_{i\circ k=0}B_k+\sum_{i\circ j=1}A_j\sum_{i\circ k=1}B_k\right)-\left(\sum_{i\circ j=0}A_j\sum_{i\circ k=1}B_k+\sum_{i\circ j=1}A_j\sum_{i\circ k=0}B_k\right) \\ &=\sum_{(j\oplus k)\circ i=0}A_jB_k-\sum_{(j\oplus k)\circ i=1}A_jB_k \\ &=FWT[C]_i \end{aligned}$

根据定义，易得 $\text{merge}$ 为 $FWT_A = \text{merge}(FWT_{A_0} + FWT_{A_1},FWT_{A_0} - FWT_{A_1})$，可以理解成左边无论右边什么最高位and 起来都是 $0$，而 $0$ 属于正贡献，所以加上右边，右边则相反。

逆变换：
$IFWT[A'] = \text{merge}(\frac{IFWT[A_0'] + IFWT[A_1']}{2}, \frac{IFWT[A_0'] - IFWT[A_1']}{2})$

逆变换的两种解释方法：根据 qyc 老师的原话，$FWT$ 的异或卷积可以看作是每一位膜 $z^2 - 1$ 的 $FFT$ 的乘法原因就是把 $0$ 看作 $0$ 次项，而 $1$ 看作一次项，那么这个多项式不断乘 $z^2$ 由于异或的性质两两抵消所以等于膜 $z^2$ 的结果为 $1$，所以就是对 $z^2 - 1$ 取模。

还有一种是根据方程组来解释，由于逆变换就是我们要把 $(a,b)$ 给还原出来 $(x,y)$，易得：$(x + y,x - y) = (a,b)$，易得 $(x,y)$ 的值，异或卷积代码：

inline void FWT_xor(int *f,int rev){for(int L = 2; L <= (1 << n); L <<= 1){for(int l = 0; l < (1 << n); l += L){for(int i = l; i < l + L / 2; ++ i){int g = f[i] * rev % mod,h = f[i + L / 2] * rev % mod;(f[i] = (g + h) % mod) %= mod;(f[i + L / 2] = (g - h + mod) % mod) %= mod; }}}
}

子集卷积

模版，求 $f(i) = \sum_{j\mid k = i,j\text{\&} k = 0} a_j b_k$

我们发现此时要求的是或是 $i$ 但是不能相交，不难发现这个第一个限制很简单就是 FWT 的板子，但是第二个有些困难，但是别忘了，我们可以在原来基础上再加一个维度表示选的子集的集合大小，具体的，设 $g(S,i) = \sum_{T\subseteq S} [popcount(T) = i] a_T$ 这个显然可以预处理，然后我们设 $f$ 是 $g$ 和 $h$ 乘起来的结果，那么我们可以让 $f(i,j) = \sum_{k = 0} ^ j g(i,k)\times h(i,j - k)$ 但是这样的满足条件的会有相交的吗，我们发现他们相交不相交并没有什么问题，但是如果他们相交，那么他们并起来的长度就不是 $i$ 了，在 IFWT 的时候自然不会被计算进去，所以我们可以这样算出来 $f$，然后再对每个可能的长度 $i\in [1,n]$ 都算一遍 IFWT，时间复杂度 $O(n ^ 22 ^ n)$

特别鸣谢 qyc 老师的教导，他解释了为什么 $|s| + |t| \rightarrow |s\mid t|$ 一定是正确的，首先我们有暴力对吧，首先做对 $A,B$ 分别做 $FWT$，这个时候一个暴力的想法是我们暴力枚举两状态 $S,T$，然后把 ${FWT_A}(S,i)\times {FWT_B}(T,j)\rightarrow {FWT_C}({S\mid T},i + j)$ 这样尽管可能 $popcount(S\mid T)\ne |S| + |T|$ 但是我们压根不会考虑这种情况，因为我们最后都是取那些 ${IFWT_C}(S,popcount(S))$ 的取计算答案，所以那些不合法的点自然会被筛掉，根据 qyc 老师的话是，我们发现最后从 ${FWT_A}(S,i)\times {FWT_B}(T,j)\rightarrow {FWT_C}({S\mid T},i + j)$ 的过程其实内部就是或卷积，所以此时直接 IFWT 是不会算重的，还有一种我自己的角度考虑，IFWT 的或卷积实际上是高维差分，所以如果不合法会在之前合法的时候被差分掉，不是很严谨感性理解即可。