多项式大乱炖

多项式大乱炖

前置知识

多项式的一些基本定义

\[f(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n \]

• 若 \(a_n\not= 0\)​，则称 \(n\)​ 是 \(f\)​​ 的次数，记作 \({\rm deg}(f)=n\)​，并称 \(a_n\)​ 为 \(f\)​ 的首项系数。

• 若 \(a_n=1\)，则称 \(f\)​ 为首一多项式。​

• 规定 \({\rm deg}(0)=-\infty\)​。

• \({\rm deg}(f(x)+g(x))\leq \max({\rm deg}(f(x)),{\rm deg}(g(x)))\)。

• \({\rm deg}(f(x)g(x))={\rm deg}(f(x))+{\rm deg}(g(x))\)。

• \([x^n]f(x)\) 表示 \(f(x)\) 中 \(x^n\)​ 项的系数。

• 多项式 \(\ne\) 多项式函数。

多项式的表示

系数表示法

对于多项式：

\[f(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n \]

称 \(\{a_0,a_1,a_2,\cdots,a_n\}\)​ 为多项式的系数表示。

点值表示法

在 \(f(x)\)​ 定义域内选择至少 \(n+1\) 个数 \(x_0,x_1,x_2,\cdots,x_n\) 代入多项式得到 \(n+1\) 个值。

称 \(\{(x_0,f(x_0)),(x_1,f(x_1)),\cdots,(x_n,f(x_n))\}\) 为多项式的点值表示。

卷积

定义

对于序列 \(a_0,a_1,\cdots,a_{n-1}\)，定义 \(\displaystyle f(x)=\sum_{i=0}^{n-1} a_ix^i\) 为这个序列的生成函数。

卷积即为生成函数的乘积在对应序列的变换上的抽象。

对于序列 \(f,g\)，其加法卷积序列 \(f\otimes g\) 满足 \(\displaystyle(f\otimes g)_k=\sum_{i=0}^kf_ig_{k-i}=\sum_{i+j=k}f_ig_j\)。

对于多项式 \(f,g\)，其多项式的卷积为 \(\displaystyle f\otimes g=\sum_{k=0}^{n+m}x^k\sum_{i+j=k}a_ib_j\)。

基本性质

• \(f\otimes g=g\otimes f\)（交换律）

通过定义不难证明。

• \((f\otimes g)\otimes h=f\otimes (g\otimes h)\)（结合律）

通过定义不难证明，两边的计算结果都为 \(\displaystyle\sum_{i+j+k=n}f_ig_jh_k\)。

• \((f\oplus g)\otimes h=(f\otimes h)\oplus (g\otimes h)\)（分配律）

其中 \((f\oplus g)_k=af_k+bg_k\)，即序列系数相加。

同样可以通过拆式子证明。

位运算卷积

定义

对于序列 \(f=\{f_0,f_1,\cdots,f_{2^n-1}\},g=\{g_0,g_1,\cdots,g_{2^n-1}\}\)，定义其位运算卷积 \(\displaystyle(f\otimes_\odot g)_k=\sum_{i\odot j=k}f_ig_j\)，其中 \(\odot\in\{\rm and,or,xor(\oplus)\}\)。

基本性质

同加法卷积，交换律，结合律，分配律均满足，证明略。

多项式乘法

定义

记 \(h(x)=f(x)g(x)\)。

则 \(\displaystyle [x^k]h(x)=\sum_{i=0}^k [x^i]f(x)[x^{k-i}]g(x)\mid k\in[0,{\rm deg}(f)+{\rm deg}(g)]\)。

本质就是多项式的加法卷积。

朴素做法

对于系数表示法而言，我们可以直接通过上述方法计算，时间复杂度：\(O(n^2)\)。

对于点值表示法而言，乘法就是把对应点点值相乘，时间复杂度 \(O(n)\)。

单位根

定义

在代数中，若 \(z^n=1|z\in \C\)，则称 \(z\) 为 \(n\)​ 次单位根。

根据复数的相关知识，不难看出 \(n\) 次单位根就是在复平面单位圆上，把单位圆 \(n\)​ 等分得到的顶点，这样的顶点一共有 \(n\) 个，且将其顺序相连可以构成一个正 \(n\) 边形。

欧拉公式：\(e^{ix}=\cos x+i\sin x\)。

根据欧拉公式，容易得到 \(n\) 个单位根为 \(e^{\frac{2\pi ki}{n}}(k\in[0,n-1])\)。

为了方便，我们设 \(\omega_n=e^{\frac{2\pi i}{n}}\)，则单位根可以表示为 \(\omega_n^i (i\in[0,n-1])\)​。

一些推广公式

• \(\omega_n^k=\cos k \dfrac{2\pi}{n}+i\sin k\dfrac{2\pi}{n}\)。
• \(\omega_{2n}^{2k}=\omega_n^k\)。
• \(\omega_n^\tfrac{n}{2}=-1\)​。
• \(\omega_n^0=\omega_n^n=1\)。

展开式子均可证明。

原根

定义与性质

对于质数 \(p\)，原根 \(g\) 满足 \(\forall i\in[0,p-1],g^i\bmod p\) 互不相同。

常见模数原根：\(998244353,1004535809,469762049\) 原根均为 \(3\)。

求原根

对 \(\varphi(p)\) 质因数分解，枚举原根判断 \(g^{\frac{\varphi(p)}{x}}\) 是否等于 \(1\) 即可。

二次剩余

定义

设 \(n\) 是正整数，若同余式 \(x^2\equiv n\pmod p\) 有解，那么称 \(n\) 是 \(\bmod p\) 的二次剩余，反之称 \(n\) 是 \(\bmod p\) 的二次非剩余。

勒让德符号

判断 \(n\) 是否为 \(p\) 的二次剩余，其中 \(p\) 为奇质数。

\[\begin{aligned} \left(\dfrac{n}{p}\right)= \begin{cases} 1&\quad n\,\texttt{是}\,p\,\texttt{的二次剩余}\\ -1&\quad n\,\texttt{是}\,p\,\texttt{的二次非剩余}\\ 0&\quad p\mid n\\ \end{cases} \end{aligned} \]

欧拉判别准则

设 \(p\) 为奇质数，对任意整数 \(n\)，有：

\[\left(\dfrac{n}{p}\right)\equiv a^{\frac{p-1}{2}}\pmod p \]

一些结论

1. \(n^2\equiv (p-n)^2\pmod p\)。

证明时展开 \((p-n)^2\) 即可。

2. 有 \(\dfrac{p-1}{2}\) 个数是 \(\bmod p\) 意义下的二次非剩余。

由结论一，可以得出 \(0\sim p-1\) 的这些数中存在 \(\dfrac{p+1}{2}\) 个不同的数，那么剩下的 \(\dfrac{p-1}{2}\) 个数就是 \(\bmod p\) 意义下的二次非剩余。

3. \((a+b)^p\equiv a^p+b^p\pmod p\)。

二项式定理展开，除了首尾两项外，其他项中的组合数都包含质因子 \(p\)，故同余式成立。

Cipolla 算法

给定 \(n,p\)，求 \(x\) 满足 \(x^2\equiv n\pmod p\)。

算法流程：

（1）根据欧拉判别准则判断原方程是否有解。

（2）随机一个数 \(a\)，使得 \(\omega^2\equiv (a^2-n)\pmod p\) 在 \(\bmod p\) 意义下是二次非剩余。

根据结论二，随机次数的期望为 \(2\) 次，复杂度有保证。

（3）找到解 \(x\equiv (a+\omega)^{\frac{p+1}{2}}\pmod p\)，由于 \(\omega^2\) 是二次非剩余，因而实数域内不存在 \(\omega\)，所以在复数域内求解即可，可以证明最后的结果合法。

证明如下：

首先证明 \(w^p\equiv -\omega\pmod p\)。

\[\omega^p=\omega\omega^{p-1}=\omega(\omega^2)^{\frac{p-1}{2}}=\omega(a^2-n)^{\frac{p-1}{2}}=-\omega \]

下面证明合法：

\[\begin{aligned} &x^2 &\\ \equiv&(a+\omega)^{p+1}&\pmod p\\ \equiv&(a+\omega)^p(a+\omega)&\pmod p\\ \equiv&(a^p+\omega^p)(a+\omega)&\pmod p\\ \equiv&(a-\omega)(a+\omega)&\pmod p\\ \equiv&a^2-\omega^2&\pmod p\\ \equiv&a^2-(a^2-n)&\pmod p\\ \equiv&n&\pmod p\\ \end{aligned} \]

牛顿迭代

主要用于求函数的零点。

对于一个函数 \(G(x)\)，求满足条件 \(G(F(x))\equiv 0\pmod{x^n}\) 的多项式 \(F(x)\)。

考虑迭代求解，假设已经求得 \(F_0(x)\) 满足：

\[G(F_0(x))\equiv 0\pmod{x^{\lceil\frac{n}{2}\rceil}} \]

将函数 \(G\) 在 \(z=F_0(x)\) 处进行泰勒展开得到：

\[G(F(x))=\sum_{i=0}^\infty\dfrac{G^{(i)}(F_0(x))}{i!}(F(x)-F_0(x))^i \]

由于 \(F(x)-F_0(x)\equiv 0\pmod{x^{\lceil\frac{n}{2}\rceil}}\)，所以 \((F(x)-F_0(x))^2\equiv 0\pmod{x^n}\)，那么把上式对 \(x^n\) 取模后只剩前两项：

\[G(F(x))\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))\equiv 0\pmod{x^n} \]

变换形式则有：

\[F(x)\equiv F_0(x)-\dfrac{G(F_0(x))}{G'(F_0(x))}\pmod{x^n} \]

根据上式求解即可。

拉格朗日插值

给定 \(n\) 个点 \((x_i,y_i)\)，确定一个 \(n-1\) 次多项式 \(f(x)\) 经过这些点。

构造如下：\(\displaystyle f(x)=\sum_{i=1}^n y_i\prod_{j\ne i}\dfrac{x-x_j}{x_i-x_j}\)。

快速傅里叶变换（FFT）

离散傅里叶变换（DFT）

把 \(n\)​ 个 \(n\)​ 次单位根代入多项式得到一个特殊的点值表示，这个过程称为离散傅里叶变换（DFT）。

离散傅里叶逆变换（IDFT）

将 DFT 得到的点值表示转化为系数表示的过程称为离散傅里叶逆变换（IDFT）。

快速傅里叶变换（FFT）

虽然 DFT 可以把多项式系数表示转化为点值表示，但时间复杂度依然是 \(O(n^2)\)​​​，考虑通过单位根的性质进行优化。

首先设多项式（注意这里为 \(n-1\)​ 次多项式，且 \(n\)​ 为 \(2\)​ 的次幂）：

\[f(x)=\sum_{i=0}^{n-1}a_ix^i=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}\mid n=2^k,k\in \N \]

按照每项次数的奇偶进行分类，再提出奇数边的 \(x\)​，不难得到：

\[f(x)=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+\cdots+a_{n-1}x^{n-2}) \]

这两边看起来非常相似，我们再设：

\[\begin{aligned} f_1(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\tfrac{n}{2}-1}\\ f_2(x)=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\tfrac{n}{2}-1}\\ \end{aligned} \]

显然有：

\[f(x)=f_1(x^2)+xf_2(x^2) \]

设 \(k<\dfrac{n}{2}\)，把 \(x=\omega_n^k\) 代入上式得到：

\[\begin{aligned} f(\omega_n^k)&=f_1((\omega_n^k)^2)+\omega_n^kf_2((\omega_n^k)^2)\\ &=f_1(\omega_n^{2k})+\omega_n^kf_2(\omega_n^{2k})\\ &=f_1(\omega_\frac{n}{2}^k)+\omega_n^kf_2(\omega_\frac{n}{2}^k)\\ \end{aligned} \]

再代入 \(x=\omega_n^{k+\tfrac{n}{2}}\)​ 得到：

\[\begin{aligned} f(\omega_n^{k+\frac{n}{2}})&=f_1((\omega_n^{k+\frac{n}{2}})^2)+\omega_n^{k+\frac{n}{2}}f_2((\omega_n^{k+\frac{n}{2}})^2)\\ &=f_1(\omega_n^{2k+n})-\omega_n^kf_2(\omega_n^{2k+n})\\ &=f_1(\omega_n^{2k})-\omega_n^kf_2(\omega_n^{2k})\\ &=f_1(\omega_\frac{n}{2}^k)-\omega_n^kf_2(\omega_\frac{n}{2}^k)\\ \end{aligned} \]

我们发现这两个多项式展开后只有符号不同，也就是说如果我们知道 \(f_1(\omega_\tfrac{n}{2}^k)\)​ 和 \(f_2(\omega_\tfrac{n}{2}^k)\)​​，就可以计算出 \(f(\omega_n^k)\) 和 \(f(\omega_n^{k+\tfrac{n}{2}})\)。

于是我们就可以递归分治 FFT 了。

时间复杂度：\(T(n)=2T(\dfrac{n}{2})+O(n)=O(n\log n)\)​。

快速傅里叶逆变换（IFFT）

结论：把原多项式的 DFT 结果作为新多项式的系数，把单位根的倒数 \(x=\{\omega_n^0,\omega_n^{-1},\cdots,\omega_n^{-(n-1)}\}\)​​ 代入新多项式，得到的每个数再除以 \(n\)，得到的就是原多项式系数。

其实就是在 FFT 的基础上取单位根的倒数再进行一次 FFT。

时间复杂度 \(O(n\log n)\)。

蝴蝶变换

为了降低递归带来的常数，我们尝试用迭代替换递归。

发现最后的序列是将原序列数字二进制水平翻转之后的结果。

所以我们只需一开始就翻转，后面用迭代模拟递归过程即可。

其中 \(\rm bit\)​​ 就是二进制位数。

三步变两步优化

设 \(A,B\) 为两个多项式，设 \(C=A+Bi\)，则有 \(C^2=A^2-B^2+2ABi\)，注意到我们要求的 \(AB\) 正好是 \(C^2\) 中虚部的 \(\dfrac{1}{2}\)，这样就只需要两次即可得到结果。

快速数论变换（NTT）

原理

FFT 利用单位根 \(\omega\) 的性质实现了分治优化多项式乘法，原根也具有这样的性质。

NTT 就是用原根代替单位根的 FFT，可以在模意义下运算，精度更高。

时间复杂度 \(O(n\log n)\)​。

多项式大杂烩

有了 NTT，我们就可以在模意义下做关于多项式的计算了。

多项式求逆

给定 \(n-1\) 次多项式 \(f(x)\)，求在 \(\bmod x^n\) 意义下的 \(g(x)\) 使得 \(f(x)\cdot g(x)\equiv1\pmod{x^n}\)。

考虑倍增，设 \(h(x)\equiv g(x)\pmod{x^{\lceil\frac{n}{2}\rceil}}\)，则有：

\[\begin{aligned} f(x)h(x)&\equiv1&\pmod{x^{\lceil\frac{n}{2}\rceil}}\\ f(x)g(x)&\equiv1&\pmod{x^{\lceil\frac{n}{2}\rceil}}\\ g(x)-h(x)&\equiv0&\pmod{x^{\lceil\frac{n}{2}\rceil}}\\ g^2(x)-2g(x)h(x)+h^2(x)&\equiv0&\pmod{x^n}\\ g(x)-2h(x)+f(x)h^2(x)&\equiv0&\pmod{x^n}\\ g(x)&\equiv h(x)(2-f(x)h(x))&\pmod{x^n}\\ \end{aligned} \]

递归算出 \(f(x)\) 在 \(\bmod x^{\lceil\frac{n}{2}\rceil}\) 意义下的逆元，按上式计算即可。

时间复杂度 \(T(n)=T(\dfrac{n}{2})+O(n\log n)=O(n\log n)\)。

多项式开根

给定 \(n-1\) 次多项式 \(f(x)\)，其中 \(x^0\) 项系数为 \(1\)，求在 \(\bmod x^n\) 意义下的 \(g(x)\)，使得 \(g^2(x)\equiv f(x)\pmod{x^n}\)。

考虑倍增，设 \(g_0^2(x)\equiv f(x)\pmod{x^{\lceil\frac{n}{2}\rceil}}\)，则有：

\[\begin{aligned} g_0^2(x)&\equiv f(x) &\pmod{x^{\lceil\frac{n}{2}\rceil}}\\ g_0^2(x)-f(x)&\equiv0 &\pmod{x^{\lceil\frac{n}{2}\rceil}}\\ (g_0^2(x)-f(x))^2&\equiv0&\pmod{x^n}\\ (g_0^2(x)+f(x))^2&\equiv4g_0^2(x)f(x)&\pmod{x^n}\\ \left(\dfrac{g_0^2(x)+f(x)}{2g_0}\right)^2&\equiv f(x)&\pmod{x^n}\\ \dfrac{g_0^2(x)+f(x)}{2g_0}&\equiv g(x)&\pmod{x^n}\\ g(x)&\equiv 2^{-1}g_0(x)+2^{-1}g_0^{-1}(x)f(x)&\pmod{x^n} \end{aligned} \]

递归算出 \(f(x)\) 在 \(\bmod x^{\lceil\frac{n}{2}\rceil}\) 意义下的平方根，求 \(g_0(x)\) 的逆，按上式计算即可。

时间复杂度 \(T(n)=T(\dfrac{n}{2})+O(n\log^2 n)=O(n\log^2 n)\)。

多项式除法

给定一个 \(n\) 次多项式 \(F(x)\) 和 \(m\) 次多项式 \(G(x)\)，求一个 \(n-m\) 次多项式 \(Q(x)\) 和一个严格小于 \(m\) 次的多项式 \(R(x)\)，使得 \(F(x)=G(x)\cdot Q(x)+R(x)\)。

设 \(F_r(x),G_r(x),\cdots\) 为 \(F(x),G(x),\cdots\) 系数翻转后的多项式，即：

\[F_r(x)=\sum_{i=0}^n a_{n-i}x^i ,\cdots \]

容易得到：\(F_r(x)=x^nF(x^{-1})\)，\(F(x)=x^nF_r(x^{-1})\)。

\[\begin{aligned} F(x)&=G(x)Q(x)+R(x)\\ F(x^{-1})&=G(x^{-1})Q(x^{-1})+R(x^{-1})\\ x^nF(x^{-1})&=x^nG(x^{-1})Q(x^{-1})+x^nR(x^{-1})\\ x^nF(x^{-1})&=x^mG(x^{-1})x^{n-m}Q(x^{-1})+x^{n-m+1}x^{m-1}R(x^{-1})\\ F_r(x)&=G_r(x)Q_r(x)+x^{n-m+1}R_r(x) \end{aligned} \]

注意到 \(R_r(x)\) 一项的系数为 \(x^{n-m+1}\)，那么把这个式子放到 \(\bmod x^{n-m+1}\) 意义下即可消除 \(R_r(x)\) 带来的影响，则：

\[\begin{aligned} F_r(x)&\equiv Q_r(x)G_r(x) &\pmod{x^{n-m+1}}\\ Q_r(x)&\equiv F_r(x)G_r^{-1}(x) &\pmod{x^{n-m+1}}\\ Q(x)&=x^nQ_r(x^{-1})&\\ R(x)&=F(x)-G(x)Q(x)&\\ \end{aligned} \]

求 \(G_r\) 的逆，按上式计算 \(Q(x),R(x)\) 即可。

时间复杂度 \(O(n\log n)\)。

多项式对数

给定一个 \(n\) 次多项式 \(f(x)\)，其中 \(x^0\) 项系数为 \(1\)，求 \(g(x)\) 满足 \(g(x)\equiv \ln f(x)\pmod{x^n}\)。

无法直接计算 \(\ln f(x)\)，考虑先求导再积分。

\[\begin{aligned} (\ln f(x))'&=\dfrac{f'(x)}{f(x)}\\ \ln f(x)&=\int\dfrac{f'(x)}{f(x)}{\rm d}x\\ g(x)&=\int\dfrac{f'(x)}{f(x)}{\rm d}x\\ \end{aligned} \]

求导，求逆元后乘起来再求积分即可。

时间复杂度 \(O(n\log n)\)。

多项式指数

给定一个 \(n-1\) 次多项式 \(g(x)\)，其中 \(x^0\) 项系数为 \(0\)，求多项式 \(f(x)\)，满足 \(f(x)\equiv e^{g(x)}\pmod{x^n}\)。

指数比较难求，考虑对两边取 \(\ln\)，变为：\(\ln f(x)\equiv g(x)\pmod{x^n}\\\)。

套用牛顿迭代得到：

\[\begin{aligned} \ln f(x)&\equiv g(x)&\pmod{x^n}\\ \ln f(x)-g(x)&\equiv 0&\pmod{x^n}\\ f(x)&\equiv f_0(x)-\dfrac{\ln f_0(x)-g(x)}{\left(\ln f_0(x)-g(x)\right)'}&\pmod{x^n}\\ f(x)&\equiv f_0(x)-\dfrac{\ln f_0(x)-g(x)}{\frac{1}{f_0(x)}}&\pmod{x^n}\\ f(x)&\equiv f_0(x)(1-\ln f_0(x)+g(x))&\pmod{x^n} \end{aligned} \]

迭代求 \(f_0(x)\)​，求对数，最后按上式乘起来即可。

时间复杂度 \(O(n\log n)\)。

多项式快速幂

给定一个 \(n\) 次多项式 \(g(x)\)，其中 \(x^0\) 项系数为 \(1\)，求多项式 \(f(x)=g^k(x)\pmod{x^n}\)。

\[f(x)=g^k(x)=e^{\ln(g^k(x))}=e^{k\ln g(x)} \]

先取 \(\ln\)​，系数乘 \(k\)​，再 \(\rm exp\) 回去即可。

时间复杂度 \(O(n\log n)\)。

多项式开根（加强版）

给定 \(n-1\) 次多项式 \(g(x)\)，不保证 \(x^0\) 项系数为 \(1\)，求在 \(\bmod x^n\) 意义下的 \(f(x)\)，使得 \(f^2(x)\equiv g(x)\pmod{x^n}\)。

就是把边界条件从赋值为 \(1\) 变为求二次剩余。

多项式幂函数（加强版）

给定一个 \(n\) 次多项式 \(g(x)\)，不保证 \(x^0\) 项系数为 \(1\)，求多项式 \(f(x)=g^k(x)\pmod{x^n}\)。

提取出来一个形如 \(ax^p\) 的公因式后把第一项的系数化为 \(1\)，再套用前面的模板即可。

快速沃尔什变换（FWT）

考虑用与 FFT 类似的方法加速位运算卷积。

求出多项式的另一种表示 \({\rm FWT}(A)\)，然后对应位置相乘，最后复原即可。

由于位运算卷积分为三种，所以分开讨论。

或卷积

设多项式 \(A,B\) 或卷积结果为 \(C\)，则有 \(\displaystyle c_k=\sum_{i\ {\rm or}\ j=k}a_ib_j\)。

构造 FWT 变换为：\(\displaystyle {\rm FWT}(A)_k=\sum_{i\ {\rm or}\ k=k}a_i\)。

下面证明 \({\rm FWT}(A)\times {\rm FWT}(B)={\rm FWT}(C)\)：

左半边：

\[\begin{aligned} &\;{\rm FWT}(A)_k\times {\rm FWT}(B)_k\\ =& \left(\sum_{i\ {\rm or}\ k=k}a_i\right) \left(\sum_{j\ {\rm or}\ k=k}b_j\right)\\ =& \sum_{i\ {\rm or}\ k=k}\sum_{j\ {\rm or}\ k=k}a_ib_j\\ =& \sum_{(i\ {\rm or}\ j)\ {\rm or}\ k=k}a_ib_j\\ \end{aligned} \]

右半边：

\[\begin{aligned} &\;{\rm FWT}(C)_k\\ =& \sum_{p\ {\rm or}\ k=k}c_p\\ =& \sum_{p\ {\rm or}\ k=k}\sum_{i\ {\rm or}\ j=p}a_ib_j\\ =& \sum_{(i\ {\rm or}\ j)\ {\rm or}\ k=k}a_ib_j\\ \end{aligned} \]

至此得证。

考虑怎样快速变换与逆变换。

注意到第一位为 \(0\) 的会对第一位为 \(1\)​ 的有贡献，\(1\) 对 \(0\) 则没有，考虑把序列分为第一位为 \(0\) 和 \(1\) 的两段递归处理后合并。

那么有 \({\rm FWT}(A)_{i+2^{n-1}}={\rm FWT}(A_0)_i+{\rm FWT}(A_1)_{i+2^{n-1}}\)。

其中 \(A_0\) 表示序列第一位为 \(0\) 的一半，\(A_1\) 表示序列第一位为 \(1\) 的一半。

而 \(A\) 的前一半就是 \(A_0\)，所以序列 \(A\) 可以通过 \(A_0\) 和 \(A_1\) 计算得到。

逆变换同理，后面一半减去前面对应的位置即可。

时间复杂度 \(O(n\log n)\)。

与卷积

设多项式 \(A,B\) 与卷积结果为 \(C\)，则有 \(\displaystyle c_k=\sum_{i\ {\rm and}\ j=k}a_ib_j\)。

类似地，构造 FWT 变换为：\(\displaystyle {\rm FWT}(A)_k=\sum_{i\ {\rm and}\ k=k}a_i\)。

此时是 \(1\) 对 \(0\) 有贡献，\(0\) 对 \(1\) 没有贡献，交换一下即可。

时间复杂度 \(O(n\log n)\)。

异或卷积

设多项式 \(A,B\) 异或卷积结果为 \(C\)，则有 \(\displaystyle c_k=\sum_{i\oplus j=k}a_ib_j\)。

异或卷积较为复杂，因为不能像前两种一样用子集关系表示。

设 \(a\odot b={\rm popcount}(a\ {\rm and}\ b)\bmod 2\)，则有：

\[(a\odot b)\oplus (a\odot c)=a\odot(b\oplus c) \]

证明时考虑当前位的取值，大力分类讨论即可。

构造 FWT 变换为：\(\displaystyle{\rm FWT}(A)_k=\sum_{i\odot k=0}a_i-\sum_{i\odot k=1}a_i\)​。

下面证明 \({\rm FWT}(A)\times {\rm FWT}(B)={\rm FWT}(C)\)：

左半边：

\[\begin{aligned} &\;{\rm FWT}(A)_k\times {\rm FWT}(B)_k\\ =& \left(\sum_{i\odot k=0}a_i-\sum_{i\odot k=1}a_i\right) \left(\sum_{j\odot k=0}b_j-\sum_{j\odot k=1}b_j\right)\\ =& \left(\sum_{i\odot k=0}a_i\sum_{j\odot k=0}b_j\right)- \left(\sum_{i\odot k=0}a_i\sum_{j\odot k=1}b_j\right)- \left(\sum_{i\odot k=1}a_i\sum_{j\odot k=0}b_j\right)+ \left(\sum_{i\odot k=1}a_i\sum_{j\odot k=1}b_j\right)\\ =& \sum_{(i\odot k)\oplus(j\odot k)=0}a_ib_j- \sum_{(i\odot k)\oplus(j\odot k)=1}a_ib_j\\ =& \sum_{(i\oplus j)\odot k=0}a_ib_j- \sum_{(i\oplus j)\odot k=1}a_ib_j\\ \end{aligned} \]

右半边：

\[\begin{aligned} &\;{\rm FWT}(C)_k\\ =& \sum_{p\odot k=0}c_p-\sum_{p\odot k=1}c_p\\ =& \sum_{p\odot k=0}\sum_{i\oplus j=p}a_ib_j- \sum_{p\odot k=1}\sum_{i\oplus j=p}a_ib_j\\ =& \sum_{(i\oplus j)\odot k=0}a_ib_j- \sum_{(i\oplus j)\odot k=1}a_ib_j\\ \end{aligned} \]

至此得证。

考虑怎样快速变换与逆变换。

同样考虑 \(k\) 和 \(k+2^{n-1}\) 的关系（以下均在 \(\bmod 2\) 意义下）。

对于 \(x<2^{n-1}\) 而言，\({\rm popcount}(x\ {\rm and}\ k)\equiv{\rm popcount}(x\ {\rm and}\ (k+2^{n-1}))\)。

对于 \(2^{n-1}\le x<2^n\) 而言，\({\rm popcount}((x-2^{n-1})\ {\rm and}\ k)\equiv{\rm popcount}(x\ {\rm and}\ k)\)，\({\rm popcount}((x-2^{n-1})\ {\rm and}\ (k+2^{n-1}))+1\equiv{\rm popcount}(x\ {\rm and}\ (k+2^{n-1}))\)。

所以有 \({\rm FWT}(A)={\rm merge}({\rm FWT}(A_1),{\rm FWT}(A_0)+{\rm FWT}(A_0)-{\rm FWT}(A_1))\)。

逆变换为 \(A={\rm merge}\left(\dfrac{A_0+A_1}{2},\dfrac{A_0-A_1}{2}\right)\)。

任意模数 NTT（MTT）

没啥用。

给定两个 \(n\) 次多项式 \(f(x),g(x)\)，求这两个多项式的乘积，系数对 \(p\) 取模。

\(n\le 10^5\)，\(a_i,b_i\le 10^9\)，\(p\le 10^9\) 且不保证能分解成 \(p=a\times 2^k+1\) 的形式。

问题的瓶颈在于：

1. 模数任意，不能直接使用 NTT。
2. 数据范围较大，不取模可能达到 \(10^{23}\)，不能直接使用 FFT。

那么相应的，我们就有两种解决方案：

1. 把 NTT 推广到任意模数的情况。
2. 通过某种办法提升 FFT 的精度。

多模数 NTT

根据第一种解决方案，我们可以找到多个符合要求的模数，在这几个模数的意义下分别做 NTT，之后再用中国剩余定理合并即可。

拆系数 FFT

把每个系数拆成 \(kM+b\) 的形式，其中 \(M\) 是常数，再分别卷积。

不难发现当 \(M=\sqrt p\) 的时候最优，此时 FFT 出来的结果是 \(10^{14}\) 级别，但需要把序列拆成四个序列，做四次 DFT，三次 IDFT，精度较低，常数较大，考虑继续优化。

优化：DFT 合并与 IDFT 合并

构造共轭式 \(P(x)=A(x)+iB(x)\)，\(Q(x)=A(x)-iB(x)\)。

推得 \(Q(\omega_n^k)\) 是 \(P(\omega_n^{n-k})\) 的共轭，那么只要 DFT 出 \(P(\omega_n^k)\)，就能得到 \(Q(\omega_n^k)\)。

然后再用 \(A(\omega_n^k)=\dfrac{P(\omega_n^k)+Q(\omega_n^k)}{2}\)，\(B(\omega_n^k)=i\dfrac{Q(\omega_n^k)-P(\omega_n^k)}{2}\)，即可通过一次 DFT 得到这两个多项式的的变换结果。

IDFT 同理，也可以把两次优化到一次，最后就变成了两次 DFT 和两次 IDFT。

多项式终极模板（vector 实现）

namespace Poly {
	const double PI = acos(-1);
	int G, IG;
	vector<int> rev, inv;
	struct comp {
		double r, i;
		comp(){}
		comp(double _r, double _i) {r = _r; i = _i;}
		comp conj() {return comp(r, -i);}
		comp operator + (const comp b) const {return comp(r + b.r, i + b.i);}
		comp operator - (const comp b) const {return comp(r - b.r, i - b.i);}
		comp operator * (const comp b) const {return comp(r * b.r - i * b.i, r * b.i + i * b.r);}
	};
	int init_root() {
		int tmp = mod - 1;
		vector<int> prm;
		for(int i = 2; i <= tmp / i; i++) {
			if(tmp % i == 0) prm.pb(i);
			while(tmp % i == 0) tmp /= i;
		}
		if(tmp > 1) prm.pb(tmp);
		for(int gg = 2; gg <= mod - 1; gg++) {
			bool flg = true;
			for(auto x : prm) {
				if(qpow(gg, (mod - 1) / x) == 1) {
					flg = false;
					break;
				}
			}
			if(flg) return gg;
		}
        assert(false); // No root
		return 0;
	}
	void poly_init(const int t) {
		G = init_root(); IG = qpow(G, mod - 2);
		inv.resize(1 << t, 0); inv[1] = 1;
		for(int i = 2; i < (1 << t); i++)
			inv[i] = inv[mod % i] * (mod - mod / i) % mod;
	}
	int extend(int n) {
		int res = 1;
		while(res < n) res <<= 1;
		return res;
	}
	int init_rev(int n) {
		int lim = extend(n);
		rev.resize(lim, 0);
		for(int i = 0; i < lim; i++) 
			rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (lim >> 1) : 0);
		return lim;
	}
	void fft(vector<comp>& A, int op, int lim) {
		A.resize(lim);
		for(int i = 0; i < lim; i++)
			if(i < rev[i]) swap(A[i], A[rev[i]]);
		for(int mid = 1; mid < lim; mid <<= 1) {
			comp wn(cos(PI / mid), op * sin(PI / mid));
			for(int len = mid << 1, pos = 0; pos < lim; pos += len) {
				comp w(1, 0);
				for(int k = 0; k < mid; k++, w = w * wn) {
					comp x = A[pos + k], y = A[pos + mid + k];
					A[pos + k] = x + y, A[pos + mid + k] = x - y;
				}
			}
		}
		if(op == 1) return;
		for(int i = 0; i < lim; i++)
			A[i].r /= lim, A[i].i /= lim;
	}
	void ntt(vector<int>& A, int op, int lim) {
		A.resize(lim);
		for(int i = 0; i < lim; i++)
			if(i < rev[i]) swap(A[i], A[rev[i]]);
		for(int mid = 2, j = 1; mid <= lim; mid <<= 1, j++) {
			int len = mid >> 1;
			int wn = qpow(G, (mod - 1) / mid);
			if(op == 0) wn = qpow(wn, mod - 2);
			for(int pos = 0, w = 1; pos < lim; pos += mid, w = 1) {
				for(int i = pos; i < pos + len; i++, w = mul(w, wn)) {
					int tmp = mul(A[i + len], w);
					A[i + len] = sub(A[i], tmp);
					A[i] = add(A[i], tmp);
				}
			}
		}
		if(op == 1) return;
		for(int i = 0; i < lim; i++)
			A[i] = mul(A[i], inv[lim]);
	}
	vector<int> poly_ntt(vector<int> A, vector<int> B) {
		int deg = sz(A) + sz(B) - 1;
		int lim = init_rev(deg);
		vector<int> C(lim);
		ntt(A, 1, lim); ntt(B, 1, lim);
		for(int i = 0; i < lim; i++)
			C[i] = mul(A[i], B[i]);
		ntt(C, 0, lim);
		C.resize(deg);
		return C;
	}
	vector<int> poly_inv(vector<int>& f, int deg) {
		if(deg == 1) return vector<int>(1, qpow(f[0], mod - 2));
		vector<int> A(f.begin(), f.begin() + deg);
		vector<int> B = poly_inv(f, (deg + 1) >> 1);
		int lim = init_rev(deg << 1);
		ntt(A, 1, lim); ntt(B, 1, lim);
		for(int i = 0; i < lim; i++)
			A[i] = mul(B[i], sub(2ll, mul(A[i], B[i])));
		ntt(A, 0, lim);
		A.resize(deg);
		return A;
	}
	vector<int> poly_dev(vector<int> f) {
		int n = sz(f);
		for(int i = 1; i < n; i++)
			f[i - 1] = mul(f[i], i);
		f.resize(n - 1);
		return f;
	}
	vector<int> poly_int(vector<int> f) {
		int n = sz(f);
		for(int i = n - 1; i; i--)
			f[i] = mul(f[i - 1], inv[i]);
		f[0] = 0;
		return f;
	}
	vector<int> poly_ln(vector<int> f, int deg) {
		vector<int> A = poly_int(poly_ntt(poly_dev(f), poly_inv(f, deg)));
		A.resize(deg);
		return A;
	}
	vector<int> poly_exp(vector<int> &f, int deg) {
		if(deg == 1) return vector<int>(1, 1);
		vector<int> B = poly_exp(f, (deg + 1) >> 1);
		B.resize(deg);
		vector<int> lnB = poly_ln(B, deg);
		for(int i = 0; i < deg; i++)
			lnB[i] = sub(f[i], lnB[i]);
		int lim = init_rev(deg << 1);
		ntt(B, 1, lim); ntt(lnB, 1, lim);
		for(int i = 0; i < lim; i++)
			B[i] = mul(B[i], add(lnB[i], 1));
		ntt(B, 0, lim);
		B.resize(deg);
		return B;
	}
	struct cp {
		int r, i;
		cp() {}
		cp(int _r, int _i) {r = _r; i = _i;}
		friend cp cpmul(cp a, cp b, int w){
			return cp(add(mul(a.r, b.r), mul(mul(w, a.i), b.i)), add(mul(a.i, b.r), mul(a.r, b.i)));
		}
	};
	int cpwr(cp a, int b, int w){
		cp res(1, 0);
		for(; b; b >>= 1, a = cpmul(a, a, w))
			(b & 1) && (res = cpmul(res, a, w), true);
		return res.r;
	}
	int cipolla(int x){
		if(qpow(x, (mod - 1) >> 1) == mod - 1) return -1;
		while(true) {
			int p = rnd(1, mod - 1);
			int w = sub(mul(p, p), x);
			if(qpow(w, (mod - 1) >> 1) == mod - 1) {
				int res = cpwr(cp(p, 1), (mod + 1) >> 1, w);
				return min(res, mod - res);
			}
		}
	}
	vector<int> poly_sqrt(vector<int>& f, int deg) {
		if(deg == 1) return vector<int>(1, cipolla(f[0]));
		vector<int> A(f.begin(), f.begin() + deg);
		vector<int> B = poly_sqrt(f, (deg + 1) >> 1);
		vector<int> IB = poly_inv(B, deg);
		int lim = init_rev(deg << 1);
		ntt(A, 1, lim); ntt(IB, 1, lim);
		for(int i = 0; i < lim; i++)
			A[i] = mul(A[i], IB[i]);
		ntt(A, 0, lim);
		for(int i = 0; i < deg; i++)
			A[i] = mul(add(A[i], B[i]), (mod + 1) >> 1);
		A.resize(deg);
		return A;
	}
	vector<int> poly_pow(vector<int> f, int k) {
		f = poly_ln(f, sz(f));
		for(auto& x : f) x = mul(x, k);
		return poly_exp(f, sz(f));
	}
	void FWTor(vector<int>& A, int op) {
		int n = sz(A);
		for(int l = 2, m = 1; l <= n; l <<= 1, m <<= 1)
			for(int j = 0; j < n; j += l) for(int i = 0; i < m; i++)
				A[i + j + m] = add(A[i + j + m], add(mul(A[i + j], op), mod));
	}
	void FWTand(vector<int>& A, int op) {
		int n = sz(A);
		for(int l = 2, m = 1; l <= n; l <<= 1, m <<= 1)
			for(int j = 0; j < n; j += l) for(int i = 0; i < m; i++)
				A[i + j] = add(A[i + j], add(mul(A[i + j + m], op), mod));
	}
	void FWTxor(vector<int>& A, bool flg){
		int n = sz(A);
		for(int l = 2, m = 1; l <= n; l <<= 1, m <<= 1){
			for(int j = 0; j < n; j += l) for(int i = 0; i < m; i++){
				int x = A[i + j], y = A[i + j + m];
				if(flg) {
					A[i + j] = add(x, y);
					A[i + j + m] = sub(x, y);
				} else {
					A[i + j] = mul(add(x, y), (mod + 1) >> 1);
					A[i + j + m] = mul(sub(x, y), (mod + 1) >> 1);
				}
			}
		}
	}
	vector<int> poly_or(vector<int> f, vector<int> g) {
		int deg = max(sz(f), sz(g));
		deg = extend(deg);
		f.resize(deg), g.resize(deg);
		FWTor(f, 1); FWTor(g, 1);
		vector<int> A(deg);
		for(int i = 0; i < deg; i++)
			A[i] = mul(f[i], g[i]);
		FWTor(A, -1);
		return A;
	}
	vector<int> poly_and(vector<int> f, vector<int> g) {
		int deg = max(sz(f), sz(g));
		deg = extend(deg);
		f.resize(deg), g.resize(deg);
		FWTand(f, 1); FWTand(g, 1);
		vector<int> A(deg);
		for(int i = 0; i < deg; i++)
			A[i] = mul(f[i], g[i]);
		FWTor(A, -1);
		return A;
	}
	vector<int> poly_xor(vector<int> f, vector<int> g) {
		int deg = max(sz(f), sz(g));
		deg = extend(deg);
		f.resize(deg), g.resize(deg);
		FWTxor(f, true); FWTxor(g, true);
		vector<int> A(deg);
		for(int i = 0; i < deg; i++)
			A[i] = mul(f[i], g[i]);
		FWTor(A, false);
		return A;
	}
} using namespace Poly;
// Remember to poly_init(Bit size)