多项式学习笔记

· · 算法·理论

前置知识：高中数学。
代码可以翻对应题目的提交记录。

多项式

相关定义

f(x)=\sum_{i=0}^na_ix_i

如上式子被我们称作关于 x 的 n 次多项式，其中 a_i 为 i 次项系数。

这个定义是我们在初中学过的，这个表示方式被我们称作“系数表示法”。But，这个我们还有别的方法可以表示这个多项式，例如图像法等，这里要介绍一种新的方法——点值表示法。

对于 n 次多项式 f，如果有 x_0,x_1\dots x_n 互不相同，那么可以由 f(x_0),f(x_1)\dots f(x_n) 唯一确定这个多项式。这个可以通过拉格朗日插值进行构造性证明。

运算

我们已知 f(x)=\sum\limits_{i=0}^na_ix_i,g(x)=\sum\limits_{i=0}^mb_ix_i(n\ge m)，那么可以定义多项式的加减法与乘除法。
加减法只需要对应系数相加减，复杂度 O(n+m) 且易证不存在复杂度更小的方法。（要是存在是不是可以得图灵奖了 hh）
至于乘除法就相对复杂。在以前的计算中，我们会把两个多项式通过乘法分配律，展开进行计算，复杂度达到了惊人的 O(nm)，且乍一看似乎也没有什么简单的优化方法。

分治乘法

我们把两个数分别按照前 \frac n2 位和 \frac n2 位进行划分，于是 (Ax+B)(Cx+D)=ACx^2+((A+B)(C+D)-AC-BD)x+BD，这玩意只需要计算 3 次乘法！于是

T(n)=3T(\frac n2)+O(n)=O(n^{\log_2^3})\approx O(n^{1.585})

（其实非常废，屁用没有）

但是，如果是在点值表示法下，乘除法会变得非常简单——直接把对应点值先相乘除嘛！于是我们产生了寻找两种表示法之间关系的想法。

我们考虑如何快速计算上面的“多点求值”&“插值”。通过前面的理论分析，我们需要得到 O(n+m) 个点值。乍一看我们只能分别求解，但是这里有一个好处：这里的点值是可以自选的。如果可以选择一些有关联的点，在求解的时候或许可以通过他们之间的联系提升效率？

DFT 帮助我们找到了有关联且方便计算的点。

DFT（离散傅里叶变换）

定义：将数组 a_0,\dots,a_{n−1}，变换得到 b_0,\dots, b_{n−1}，设 \omega_n是 n 次单位根

b_j=\sum_{i=0}^{n-1}a_i\omega_n^{ij}

事实上，这玩意的原始用途并不是多项式，而是用于信号处理。不过好巧不巧，这玩意在多项式上有很大的作用。如果我们把 a 看作多项式，则 b_j=a(\omega_n^j)，即上面所说的“点值”。

但是 DFT 本身并没有加快程序的速度，仍然是 O(n^2) 的算法。

IDFT（离散傅里叶逆变换）

定义：将数组 b_0,\dots,b_{n−1}，变换得到 a_0,\dots, a_{n−1}，设 \omega_n是 n 次单位根

a_i=\frac1n\sum_{j=0}^{n-1}b_j\omega_n^{-ij}

证明？带入即可。

单位根反演

由 IDFT 发现的一个重要结论：

[n\ |\ k]=\frac1n\sum_{i=0}^{n-1}\omega_n^{ik}

证明如下：

这个东西可以方便判断倍数，例如 P5591 小猪佩奇学数学。

FFT（快速傅里叶变换）

当 n=2^k 时，上面的 DFT 其实可以分治处理，复杂度变为 O(n\log n)，即为 FFT。

设 n=2m，首先我们对多项式 A(x) 进行拆分：

A_0(x)=a_0+a_2x+a_4x^2\dots a_{n-2}x^{m-1} A_1(x)=a_1+a_3x+a_5x^2\dots a_{n-1}x^{m-1} A(x)=A_0(x^2)+x\cdot A_1(x^2)

其实就是把奇偶分开了。
我们观察这样做有什么好处。设 i\in[0,m)，可以得出结论：

\begin{aligned} A(\omega_n^i)&=A_0(\omega_n^{2i})+\omega_n^i\cdot A_1(\omega_n^{2i})\\ &=A_0(\omega_m^i)+\omega_n^i\cdot A_1(\omega_m^i)\\ A(\omega_n^{i+m})&=A(-\omega_n^i)\\ &=A_0(\omega_n^{2i})-\omega_n^i\cdot A_1(\omega_n^{2i})\\ &=A_0(\omega_m^i)-\omega_n^i\cdot A_1(\omega_m^i) \end{aligned}

我们发现，通过将原问题划分为两个子问题，可以在线性复杂度内找到 A(x) 在每个 n 次单位根位置的值。即 T(n)=2(\dfrac n2)+O(n)=O(n\log n)。

解决了求值，别忘了还有插值。显然插值对应了求值的逆运算，所以可以用 IDFT 解决。

Code

注意，在理想状态下，算完之后的虚部应当为 0。不过既然是小数运算，必然有误差，所以应当对实部进行四舍五入。

void FFT(Complex a[], int tot, int step, int inv)//inv表示正/逆变换 
{
    if(n == 1) return;
    int m = n >> 1;
    FFT(a, m, step << 1, inv, a);
    FFT(a + step, m, step << 1, inv, b + m);
    Complex w = (Complex){cos(PI / m), sin(PI / m) * inv}, wk = (Complex){1, 0};

    for(int i = 0; i < m; i++, wk = wk * w)
    {
        Complex x = a[i], y = a[i + m] * wk;
        a[i] = x + y, b[i + m] = x - y;
    }
}

int main()
{
    FFT(a, tot, 1, 1), FFT(b, tot, 1, 1);
    for(int i = 0; i < tot; i++) a[i] = a[i] * b[i];
    FFT(a, tot, 1, -1);
    for(int i = 0; i < tot; i++) a[i].x /= n, a[i].y /= n;
    return 0;
}

然而这玩意根本过不了板题

上述代码的常数很大，我们尝试进行优化。
观察 FFT 的前五行都在干些什么？其实是以一个奇怪的顺序进行 swap 操作。稍加思考可以发现，最后就是把 a_i 变为了 a_{rev(i)}，其中 rev 表示二进制下的回文串。

而且 \omega 的计算也很冗余，这里把同一个 \omega 计算了很多次，这不是我们期望的，于是我们决定用数组将单位根存储下来。此外，double 的除法相对于乘法而言超级慢，所以程序后面 \div n 操作非常不友好，换成乘法即可。

事实上，优化完的板子由于运算次数减少了，不仅常数会减少，精度也会变高！

常数优化

常见的 fft 有两个比较通用且显著的优化：

三次变两次

上面的程序中，每进行一次多项式乘法都会执行 3 次 fft。有没有什么减少 fft 执行次数的方式呢？
注意到两个多项式在赋初值的时候，他们的虚部都是空着的，这造成了浪费，于是我们考虑把 B 这个多项式放到 A 的虚部上。于是原先的 A\times B 就变成了 A^2=(a+bi)^2=a^2-b^2+2abi。我们惊奇的发现，虚部的一半刚好就是我们要的答案！（可惜的是这个优化不能放到 ntt 上用了）

长度优化

由于最后的长度必然是 2^k，所以只要优化了长度，那么至少都可以省去 \frac12 的常数。所以说，平时在进行多项式乘法的时候，最好是精打细算一下长度。
当然这里并不是仅仅指避免开多余的长度——就算是不多于的长度，我们也可以压缩。举个例子，(1+x+2x^2)(-1+x-x^2)=-1-2x^2+x^3-2x^4，但正常而言需要开 tot=8 项，比较浪费。如果把 tot\gets4，会发生什么呢？结果是：-3-2x^2+x^3。

上面的例子并不是巧合，这是有迹可循的：4\sim7 这些项被我们扔掉了，那么他们就会叠加到 0\sim3 项上。更普通的：如果结果的长度应该为 n 项，而只开了 \frac n2 项，导致的结果就是后面一半的系数会叠加到前面一半身上！（读者自证不难）
因此，如果我们用不到那些被影响的项，完全可以压缩长度。举个例子，我们把 2n 压缩成了 n，而答案的上限是 \frac32n 项，那么我们这样做之后，\frac12n\sim n 这些项我们依旧可以正常访问！
（比较 nice 的是，这个优化显然可以丢到 ntt 上）

NTT（快速数论变换）

FFT 由于涉及到了小数运算，带来了两个缺陷：其一是常数问题，不过也还能接受，毕竟慢得不会太多；其二是精度问题，这就很头疼了，因为避免不了！于是，在涉及到“取模”的问题时，FFT 就显得力不从心了。

带来小数运算的根本原因是单位根！那么有没有什么东西，可以代替掉单位根呢？
由于 \{\omega_n^0,\omega_n^1\dots\omega_n^{n-1}\} 构成了一个乘法群，我们从群论的视角思考单位根的性质：

似乎单位根就用到了这些性质？
那有没有替代品呢？

在一般情况下其实不太好找，不过当模数为质数的时候，有一个非常妙的东西——原根！容易发现，当模数为质数的时候，原根满足上述单位根的所有性质！换句话说，我们可以用原根来代替单位根。

NTT 模数

在上面代码中用到了这样一个单位根：(\cos\frac\pi m, \sin\frac\pi m)=(\cos\frac{2\pi}{2m}, \sin\frac{2\pi}{2m})，如果翻译为原根，会变成 g^{\frac{mod-1}{2m}}，其中 g 为原根。于是我们发现，这里要求 2m|mod-1，即 mod=t\cdot2^k+1 且 k 要足够大！（例如要求 k\ge20）

这便是 NTT 模数的由来。因为经过上面的推导，我们发现并非所有数都可以做 ntt。

常见的 NTT 模数：

\begin{aligned} 998244353&=119\times2^{23}+1&g=3\\ 1004535809&=479\times2^{21}+1&g=3\\ 950009857&=453\times2^{21}+1&g=7\\ 1811939329&=27\times2^{26}+1&g=13\\ 469762049&=7\times2^{26}+1&g=3\\ 167772161&=5\times2^{25}+1&g=3 \end{aligned}

一个有趣的故事

事实上，在以前的 OI 中，出题人偏爱的模数只有 10^9+7 一个。后来出现了 NTT，题目中需要有 NTT 模数，于是出现了一部分题目的模数为 998244353。于是出现了这样的情况：只要看到 998244353，这个题必然是 NTT！
为了防止这种情况发生，有人便开始呼吁：我们以后都把题目的模数改为 998244353，就不会发生这种情况了。也正因如此，导致了 10^9+7 和 998244353 两个数遍布题目的现状。

小技巧

由于取模运算很慢，且足足执行了 O(n\log n) 次，我们考虑优化。在下面的模板中采用短路取模进行优化，获得了不错的效果。不过，我们能不能从根本上解决问题，直接不取模呢？事实上，答案是可以的。

void NTT(int a[], int inv)
{
    if(!rev[bit][1])
        for(int i = 0; i < tot; i++) rev[bit][i] = rev[bit][i >> 1] >> 1 | (i & 1) << bit - 1;
    static unsigned LL t[N];
    for(int i = 0; i < tot; i++) t[i] = a[rev[bit][i]];
    static int g[N];
    g[0] = 1;
    for(int len = 1; len < tot; len <<= 1)
    {
        int g1 = ksm(~inv ? 3 : 332748118, (mod - 1) / (len << 1));
        for(int i = len - 2; i >= 0; i -= 2)
            g[i] = g[i >> 1], g[i + 1] = 1ll * g[i] * g1 % mod;
        for(int i = 0; i < tot; i += len << 1)
            for(int j = i; j < i + len; j++)
            {
                int g1 = t[j + len] * g[j - i] % mod;
                t[j + len] = t[j] - g1 + mod, t[j] += g1;
            }
    }
    for(int i = 0; i < tot; i++) a[i] = t[i] % mod;
    if(inv == -1)
    {
        int t = ksm(tot, mod - 2);
        for(int i = 0; i < tot; i++) a[i] = 1ll * a[i] * t % mod;
    }
}

在这里我们采取了 unsigned longlong 进行处理，且相比于下面的模板有微小的加速。事实上，我也用普通的 longlong 处理过这个问题，不过很遗憾 WA 飞了。问题在哪里呢？
我们考虑值域变化：注意到代码中的 t 数组，一开始是 [0,mod) 的，在进行完单次变换后，大小来到了 [0,2mod)，再次变换会变为 [0,3mod)，即每次值域 +mod。由于总共会变换 bit 次，因而我们的值域会达到 [0,(bit+1)\cdot mod)。其实在程序中最令人头疼的是 t\cdot g 的那步运算，t 的上限为 bit\cdot mod，故程序的运算上界为 bit\cdot mod^2。 当 mod=998244353 时，longlong 的范围为 9\cdot mod^2，显然不够用，而 unsigned longlong 为 18\cdot mod^2，可以解决 bit\le18，即 tot\le262144 时的问题。

当然，上述分析仅仅是理论上限，真正想达到这个上界不大可能，所以这个模板平时不会炸，但依然不建议平时用，也就卡常的时候掏出来用用得了……
（文章末尾有一些更高端的卡常技巧，这个其实真的没啥用）

第二类斯特林数·行

学会了 NTT，就可以解决这个经典问题了。

事实上，第二类斯特林数存在通项公式：

S(n,m)=\sum_{i=1}^m\frac{(-1)^{m-i}i^n}{i!(m-i)!}

证明如下：

考虑如下问题：
有一个长度为 n 的序列且满足 a_i\in[1,m]，要求 1\sim m 的每个数至少出现一次，求方案数？

【法一】原问题等价于先将序列分为 m 个部分，再依次标号 1\sim m，所以方案数为 S(n,m)\cdot m!

【法二】考虑二项式反演。既然要求每个数都要出现，容斥即可，方案数为 \sum_{i=1}^m(-1)^{m-i}\cdot C_m^i\cdot i^n。

于是联立上面的式子，将组合数用阶乘表示，即可得证。

我们将上面的式子整理可得：

S(n,m)=\sum_{i=1}^m\frac{(-1)^{m-i}}{(m-i)!}\cdot\frac{i^n}{i!}

于是构造 a_i=\dfrac{(-1)^i}{i!},b_i=\dfrac{i^n}{i!}，就可以使用多项式乘法求解了。

分治 NTT

一个小技巧。

先讲讲比较简单的版本。计算:

\prod_{i=1}^n(x-a_i)

如果依次相乘，每次乘法的规模是 O(n) 的，所以复杂度为 O(n^2\log n)。问题在于 NTT 的复杂度与长度相关，于是我们考虑优化多项式相乘的长度，直接分治即可。这样做的复杂度会被优化到 O(n\log^2n)。

但是如果是上面的模板题那样，前后具有依赖关系，好像就不能朴素分治了，需要改一下。其实只需要用类似 CDQ 分治的思路即可，即在之前的基础上统计前面一半区间对后一半的贡献。

void cdq(int l, int r)
{
    if(l == r) return;
    int mid = l + r >> 1;
    cdq(l, mid);

    bit = lg2(r - l - 1 + mid - l) + 1, tot = 1 << bit;
    for(int i = 0; i < tot; i++) a[i] = b[i] = 0;
    for(int i = l; i <= mid; i++) a[i - l] = f[i];
    for(int i = 1; i <= r - l; i++) b[i - 1] = g[i];
    NTT(a, 1), NTT(b, 1);
    for(int i = 0; i < tot; i++) a[i] = 1ll * a[i] * b[i] % mod;
    NTT(a, -1);
    for(int i = mid + 1; i <= r; i++) f[i] = (f[i] + a[i - l - 1]) % mod;

    cdq(mid + 1, r);
}

总用时 1.05s，还不错，不过用下面的多项式求逆解决这道题只需要 250ms。

具体地说，显然题目告诉我们存在关系 f=fg。不过直接这么写是错的，因为显然这个方程等价于 (g-1)f=0，没有什么意义。问题在于初值 f_0=1，我们没有体现到方程上，于是方程应该为 f=fg+1，解得 f=\dfrac1{(1-g)}。所以只要我们会求逆，这道题就迎刃而解了。而求逆的复杂度为 O(n\log n)，所以快上不少。

多项式全家桶

多项式求逆

首先要明白一件事：任意一个有限次数多项式，它的逆元一定是无限次数的！这一点很容易说明：采用反证，如果 n 次多项式的逆元是 m 次的，那么他们的乘积一定是 n+m 次的多项式而非常数 1。
所以说我们看到题目中的“\pmod{x^n}”，其实是让你求逆元的前 n 项而已。

考虑倍增。在后面的很多推导中都会这样。

\begin{aligned} G&\equiv\frac1F&\pmod{x^n}\\ G'&\equiv\frac1F&\pmod{x^{2n}}\\ G'-G&\equiv0&\pmod{x^n}\\ G'^2-2G'G+G^2&\equiv0&\pmod{x^{2n}}\\ G'-2G+G^2F&\equiv0&\pmod{x^{2n}}\\ G'&\equiv2G-G^2F&\pmod{x^{2n}}\\ \end{aligned}

倍增有一个很重要的好处：虽然倍增自带 \log，每一层的复杂度也带 \log，但总复杂度仍然为 O(n\log n)，这一点可以依靠主定理求得。

多项式开方

注意到题目中含有条件 a_0=1。这个限制其实是为了避免二次剩余，考虑极限情况：这个多项式的次数为 0，那么就是在求二次剩余。所以多项式开根（加强版）只需要多一个二次剩余的板子即可。

同样倍增。

\begin{aligned} G&\equiv\sqrt F&\pmod{x^n}\\ G'&\equiv\sqrt F&\pmod{x^{2n}}\\ G'-G&\equiv0&\pmod{x^n}\\ G'^2-2G'G+G^2&\equiv0&\pmod{x^{2n}}\\ F-2G'G+G^2&\equiv0&\pmod{x^{2n}}\\ G'&\equiv\frac{F+G^2}{2G}&\pmod{x^{2n}}\\ \end{aligned}

多项式 ln

首先注意到题目中的限制 a_0=1，事实上，当且仅当满足这个限制时，多项式在模意义下存在对数函数。也同样的，当且仅当 a_0=0 时，存在指数函数。
证明

明确了这一点之后，其实对数函数的推导本身不算太难。

\ln F=\int(\ln F)'=\int\frac{F'}F

需要注意的是，因为是不定积分，所以最后算出来的常数项应该为 0，写代码的时候要注意一下。

多项式 exp

牛顿迭代

考虑这样一个问题：已知一个可导单调函数 f(x)，求解函数的零点。
或许会想到二分。的确，二分可以在 \log 复杂度内解决这个问题，不过这里将引入一个更加迅速的方法：牛顿迭代。

考虑当前答案为 t。由于是可导函数，我们可以画出 f 在 t 处的切线：y=f'(t)x+f(t)-f'(t)t。接着，我们令 t' 为这个切线方程的零点，即 t'=t-\dfrac{f(t)}{f'(t)}。容易想到，t' 一定在逼近函数的零点，事实上，依靠泰勒展开可以证明，这样迭代一次可以使精度翻倍。显然，这比起二分快得多。

依靠牛顿迭代求解 exp

\begin{aligned} G&\equiv e^F&\pmod{x^n}\\ \ln G&\equiv F&\pmod{x^n}\\ \ln G-F&\equiv 0&\pmod{x^n}\\ \end{aligned}

即：我们要找到一个多项式 G，使得 G 为函数 \ln x-F 的零点。

同样倍增（倍增的是精度），因为可以使用牛顿迭代了。(在套用牛顿迭代的时候，我们 x 是变量，所以这里的多项式也要视为变量，即求导后为 1)

\begin{aligned} \ln G-F&\equiv 0&\pmod{x^n}\\ H&\equiv G-\frac{\ln G-F}{(\ln G-F)'}&\pmod{x^{2n}}\\ H&\equiv G-\frac{\ln G-F}{\dfrac1G}&\pmod{x^{2n}}\\ H&\equiv G(1-\ln G+F)&\pmod{x^{2n}}\\ \end{aligned}

多项式快速幂

显然当 k 比较小的时候可以朴素快速幂，复杂度 O(n\log n\log k)。
然后发现题目保证的是 k 的位数。
事实上，对 k 使用扩展欧拉定理取模之后依旧可以朴素快速幂，不过带 \log，我们不喜欢。

显然有 A^k=\exp(\ln(A^k))=\exp(k\ln(A))，于是……
复杂度 O(n\log n)。

bbbbbut，难道就完了吗？ 打开加强版一看，没有了 a_0=1 的限制。注意到前面的 \ln 要求 a_0=1，\exp 要求 a_0=0。（满足 \ln 的就可以了，后面 \exp 显然满足）
其实方法很简单，把首位强制变为 1 即可，即让整个多项式 \div a_0。不过又出问题了，a_0=0 会寄飞，所以还要先提取 0。
最后注意一下，由于最后要把 a_0^k 乘回去，所以这里的 k 不同于上面，是对 \varphi(mod) 取模的。

多项式除法/取模

有点困难

这里要用到一点人类智慧（其实数学上还是非常套路的）：对于多项式 f(x)=\sum_{i=0}^n a_ix^i，我们构造 f'(x)=\sum_{i=0}^na_{n-i}x^i，即将系数进行 reverse 操作，那么有 f'(x)=x^nf(\frac1x)。

由题，我们有 F(x)=Q(x)G(x)+R(x)，其中各多项式的次数分别为 n,m,n-m,m-1。于是开始施法：（和上面一样，我们用 f' 表示系数的 reverse 操作）

\begin{aligned} F(x)&=Q(x)G(x)+R(x)\\ x^nF(x)&=x^mQ(x)x^{n-m}G(x)+x^{m-1}R(x)\cdot x^{n-m+1}\\ F(x)'&=Q(x)'G(x)'+R(x)'\cdot x^{n-m+1}\\ F(x)'&\equiv Q(x)'G(x)'\pmod{x^{n-m+1}}\\ \end{aligned}

我们可以根据这个式子求出 G(x)'\pmod{x^{n-m+1}}。注意，由于 G(x) 本来就是 n-m 次多项式，所以取模其实根本没有影响。因此，我们用 O(n\log n) 求出了 G(x)。再带回原式，R(x)=F(x)-Q(x)G(x) 即可求得。

模板中的 a,r 可以相同。

常系数齐次线性递推

喵呜\~\~\~来咯\~\~\~

为了后续方便，我们令首项下标为 0。

以前的做法是矩阵快速幂，复杂度 O(k^3\log n)，一般适用于 k 比较小的情况，最经典的是斐波那契。不过遇到这种情况，矩阵当然是行不通了，考虑使用多项式代替矩阵。

我们都知道，对于 k 阶线性递推，存在 k 次的特征方程。我们从这个特征方程下手，具体的，我们令 f(0)=1,f(1)=x,f(t)=x^t。我们考虑用 x^0\sim x^{k-1} 表示其他值，显然根据特征方程，这点一定可以达到。
于是操作如下：

f(t)=x^t=x\cdot x^{t-1}=\sum_{i=1}^kp_ix^i=p^kx^k+\sum_{i=1}^{k-1}p_ix^i

于是我们成功通过递推的方式表示出了 f(t)。显然，单次递推的复杂度为 O(k)，所以复杂度为 O(nk)。（和暴力一样）

不过，以上做法给了我们灵感：既然 f(t)=x^t，那么我们能否使用快速幂的方式计算呢？显然是可以的。
于是，我们将上面的操作总结为：先把两个 k-1 次多项式相乘，再将结果降次至 k-1 次。

这样做的话，暴力乘法和降次的复杂度均为 O(k^2)，故复杂度 O(k^2\log n)，比矩阵快速幂优秀不少。

不过还不够，这道题的数据范围并不支持我们使用 k^2 的计算方式。考虑优化，其实非常简单，多项式乘法可以使用 NTT 优化，而降幂的本质其实就是对特征方程取模，可以用多项式取模优化。于是复杂度就变为了 O(k\log k\log n)。

另外，如果使用快速幂会带有一个 2 的常数，但如果采用倍增可以避免。

f[n] = 1;
for(int i = 0; i < n; i++) f[i] = (mod - b[n - i]) % mod;

res[0] = 1;
for(int i = 32; ~i; i--)
{
    bit = lg2(n - 1 << 1) + 1, tot = 1 << bit;
    for(int i = n; i < tot; i++) res[i] = 0;
    NTT(res, 1);
    for(int i = 0; i < tot; i++) res[i] = 1ll * res[i] * res[i] % mod;
    NTT(res, -1);
    Div(res, tot - 1, f, n, t, res);
    if(m >> i & 1)
    {
        for(int j = n; j; j--) res[j] = res[j - 1];
        res[0] = 0;
        for(int j = 0; j < n; j++) res[j] = (res[j] + 1ll * res[n] * b[n - j]) % mod;
    }
}
for(int i = 0; i < n; i++) ans = (ans + 1ll * res[i] * a[i]) % mod;

多项式多点求值

有一个范德蒙德转置矩阵的做法，是集训队论文中提到的，常数比这里的更小，但不做讲解。

考虑 n 很大，但 m 很小怎么做？
由于暴力带入的复杂度达到了 O(nm)，我们考虑降幂。注意到 g(x)=\prod(x-a_i) 对于所有 x=a_i 都为 0，所以我们可以让原多项式 f(x) 对 g(x) 取模，再带入 a_i，易证答案不变。这种方法的复杂度为 O(n\log n+m^2)。
但是 m 也很大怎么办呢？注意到我们完全没有必要全部一起求，比如可以将 a_i 分块，一块一块求，如果块长为 \sqrt m，则复杂度为 O(\sqrt m(n\log n+m))，好像还不错。

事实上，我们可以采用分治的方法。令 g(l,r)=\prod\limits_{i=l}^r(x-a_i)，假设当前的分治区间为 [l,r]，多项式为 f(l,r)，于是可以得到 f(l,mid)=f(l,r)\bmod g(l,mid),f(mid+1,r)=f(l,r)\bmod g(mid+1,r)，然后继续分治即可。当 l=r 时，显然 f(l,r) 为 0 次，它的常数项即为 f(a_l) 的值。
这样做的复杂度为 O(n\log n)。不过，虽然多项式除法的复杂度为 O(n\log n)，但是其实是有下界的。举个例子，我们会调用多项式求逆，而在求逆的时候会求常数项的逆元，因此下限至少是 O(n\log n+\log mod)。而在 n 很小的时候，可能后面那些东西的复杂度更高了，导致代码常数其实很大。建议在 r-l\le150 时直接暴力，实测时间为原来的 \frac12。

在实现的时候，可以用分治 NTT 求得上面的 g 数组。

void init(int l, int r, int deep)
{
    if(l == r)
    {
        g[deep][l] = mod - a[l];
        return;
    }
    int mid = l + r >> 1;
    init(l, mid, deep + 1), init(mid + 1, r, deep + 1);
    static int x[N], y[N];
    for(int i = l; i <= mid; i++) x[i - l] = g[deep + 1][i];
    for(int i = mid + 1; i <= r; i++) y[i - mid - 1] = g[deep + 1][i];
    x[mid - l + 1] = y[r - mid] = 1;
    bit = lg2(r - l + 1) + 1, tot = 1 << bit;
    for(int i = mid - l + 2; i < tot; i++) x[i] = 0;
    for(int i = r - mid + 1; i < tot; i++) y[i] = 0;
    NTT(x, 1), NTT(y, 1);
    for(int i = 0; i < tot; i++) x[i] = 1ll * x[i] * y[i] % mod;
    NTT(x, -1);
    for(int i = l; i <= r; i++) g[deep][i] = x[i - l];
}

void solve(int l, int r, int deep)
{
    if(r - l <= 150)
    {
        int m = r - l;
        for(int i = l; i <= r; i ++ )
        {
            b[i] = 0;
            for(int j = m; j >= 0; j -- )
                b[i] = (1ll * b[i] * a[i] + f[deep][j]) % mod;
        }
        return;
    }
    int mid = l + r >> 1;
    static int x[N];
    for(int i = l; i <= mid; i++) x[i - l] = g[deep + 1][i];
    x[mid - l + 1] = 1;
    Div(f[deep], r - l, x, mid - l + 1, f[0], f[deep + 1]);
    solve(l, mid, deep + 1);
    for(int i = mid + 1; i <= r; i++) x[i - mid - 1] = g[deep + 1][i];
    x[r - mid] = 1;
    Div(f[deep], r - l, x, r - mid, f[0], f[deep + 1]);
    solve(mid + 1, r, deep + 1);
}

拉格朗日插值

前置知识：拉格朗日插值。
谁规定了自己的前置知识不能是自己

带修

实际上就是重心拉格朗日的应用。

f(x)=\sum_iy_i\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}

这个是朴素的插值。

f(x)=\prod_i(x-x_i)\cdot\sum_i\frac{y_i}{(x-x_i)\prod_{j\ne i}(x_i-x_j)}

这是重心拉插。

于是我们只需要维护 \prod_{j\ne i}(x_i-x_j)，每次查询时计算 \prod_i(x-x_i) 即可。需要注意，如果 x 与某个 x_i 相等，那么程序就会除 0，所以需要特判。
代码很简单，懒得贴了。

多项式快速插值

在上面的公式中，既然我们可以把分子提出来，为什么不能把分母提出来呢？我们试试转化分母：

x\in x_i,\prod_{j\ne i}(x-x_j)=\frac{\prod\limits_jx-x_j}{x-x_i}

后面那坨分子分母都为 0，于是洛必达：

\prod_{j\ne i}(x-x_j)=\left.\left(\prod\limits_jx-x_j\right)'\right|\large_{x=x_i}

（服了，这 \LaTeX 怎么这么丑……）

不过相比之前的东西起码好算得多了。我们令 g(x)=\prod\limits_i(x-x_i)，那么根据之前的结论：

f(x)=\sum_i\frac{y_i}{g'(x_i)}\prod\limits_{j\ne i}(x-x_j)

分母好算多了，直接套用多点求值即可。

剩余的部分可以分治 NTT 了计算了，好耶！
复杂度同多点求值，O(n\log^2n)。

多项式复合逆

来点黑科技。该部分内容适合多项式老油条食用，新手建议跳过。
前置知识：拉格朗日反演，Bostan-Mori。

我们知道 F(G(x))=x，于是有 F^k(G(x))=x^k。根据扩展拉格朗日反演，有结论：

[x^n]F^k=\frac1n[x^{n-1}]kx^{k-1}\left(\dfrac xG\right)^n

大力推式子，尝试用 F^k 去表示 G。

\frac nk[x^n]F^k=[x^{n-k}]\left(\dfrac xG\right)^n\\ \sum_{k=1}^nx^{n-k}\frac nk[x^n]F^k=\left(\dfrac Gx\right)^{-n}\\ G=x\left(\sum_{k=1}^nx^{n-k}\frac nk[x^n]F^k\right)^{-\frac1n}

设 H=\sum\limits_{k=1}^nx^{n-k}\frac nk[x^n]F^k，那么 H^{-\frac1n} 可以多项式快速幂解决。有一个小细节需要注意的是，在算快速幂的过程中需要求常数项的 n 次剩余。这个问题本身比较困难，但在这里 h_0=[x^n]F^n=f_1^n，于是 n 次剩余直接就是 f_1。为了方便，可以直接把常数项提出来。具体的，设 v=f_1，那么上述式子会变形为：

G=\frac xv\left(\sum_{k=1}^nx^{n-k}\frac n{kv^{n-k}}[x^n]\left(\dfrac Fv\right)^k\right)^{-\frac1n}

剩下的难点是如何计算所有的 [x^n]F^k。用二元生成函数的形式去刻画，那么我们需要求解 [x^n]\frac1{1-yF}。其中，保证 f_0=0。
将 x 视作主元，y 视作系数，那么这个问题本质就是在提取分式型形式幂级数的远端系数。因此， Bostan-Mori 可以在 O(n\log^2n) 时间复杂度求解。这里可能需要理解一下，为什么此处二元函数跑 Bostan-Mori 的复杂度，与朴素的一元函数相同。对于第 k 层而言，我们只需要保留 F(x) 的 O(\frac n{2^k}) 项，而每一项的系数是一个 O(2^k) 次关于 y 的多项式。因此每一层有 O(n) 个值，总共 \log n 层，每层跑多项式乘法还有一个 \log，总复杂度为 O(n\log^2n)。
最后补充一下二元多项式的乘法。我们用 F(n,m) 表示 \sum\limits_{i=0}^n\sum\limits_{j=0}^mf_{i,j}x^iy^j，假设我们要处理 F(a,b)\times G(c,d)。考虑将其映射为一个一元多项式，f_{i,j}\to[z^{i(b+d+1)+j}]，这是因为 y 的最高次数只有 b+d，所以这样是对的。

至此，我们在 O(n\log^2n) 复杂度内解决了多项式复合逆问题。

卡常宝藏！

气温大哥稳定发力。

还有指令集大手子写的 多项式超级模板。

FWT（快速沃尔什变换）

本质上不是多项式，所以没有丢进上面的封装里。

前置知识：高维前缀和。

FWT_or

\begin{aligned} C_i&=\sum_{j\cup k=i}A_jB_k\\ \sum_{j\subseteq i}C_j&=\sum_{j\cup k\subseteq i}A_jB_k\\ \sum_{j\subseteq i}C_j&=\sum_{j\subseteq i}A_j\sum_{j\subseteq i}B_j \end{aligned}

然后就可以做了。发现这个东西其实就是高维前缀和，于是 IFWT_or 本质上就是高维差分。

void FWT_or(int a[], int n, int inv)
{
    for(int i = 1; i < n; i <<= 1)
        for(int j = 0; j < n; j += i << 1)
            for(int k = j; k < j + i; k++)
                ~inv ? ((a[k + i] += a[k]) >= mod && (a[k + i] -= mod)) : ((a[k + i] -= a[k]) < 0 && (a[k + i] += mod));
}

FWT_and

\begin{aligned} C_i&=\sum_{j\cap k=i}A_jB_k\\ \sum_{i\subseteq j}C_j&=\sum_{i\subseteq j\cup k}A_jB_k\\ \sum_{i\subseteq j}C_j&=\sum_{i\subseteq j}A_j\sum_{i\subseteq j}B_j \end{aligned}

与运算和或运算挺像的，所以把前缀换成后缀就可以解决。

void FWT_and(int a[], int n, int inv)
{
    for(int i = 1; i < n; i <<= 1)
        for(int j = 0; j < n; j += i << 1)
            for(int k = j; k < j + i; k++)
                ~inv ? ((a[k] += a[k + i]) >= mod && (a[k] -= mod)) : ((a[k] -= a[k + i]) < 0 && (a[k] += mod));
}

FWT_xor

前两个只是开胃菜。这个怎么解决呢？
（尝试发明高维中缀和ing……）

我们观察前面两个的代码——好像和 NTT 极其相似，唯一的区别是括号里面的东西。我们发现，前面的 \operatorname{or} 和 \operatorname{and} 其实还是在找一个“点值”，然后直接乘除。那么这里的 \operatorname{xor} 怎么寻找点值呢？

以 n=2 为例：

我们发现无论是前缀和还是后缀和，都是为了因式分解，C_0+C_1=(A_0+A_1)(B_0+B_1)。不过光这一个式子只能还原出 FWT_or 和 FWT_and，无法得到 FWT_xor，所以还需要一个 C_0-C_1=(A_0-A_1)(B_0-B_1)

此时就得到点值了，即把 (a,b) 变换为 (a+b,a-b)。不过在 IFWT_xor 的时候，我们发现需要 \div2，所以我们把所有的除法提出来，总共就是 \div n。（每一维都要 \div2）

void FWT_xor(int a[], int n, int inv)
{
    for(int i = 1; i < n; i <<= 1)
        for(int j = 0; j < n; j += i << 1)
            for(int k = j; k < j + i; k++)
            {
                int x = a[k], y = a[k + i];
                (a[k] = x + y) >= mod && (a[k] -= mod), (a[k + i] = x - y) < 0 && (a[k + i] += mod);
            }
    if(inv == -1)
    {
        int t = ksm(n, mod - 2);
        for(int i = 0; i < n; i++) a[i] = 1ull * a[i] * t % mod;
    }
}

小结

我们发现这些代码与 NTT 都极其相似，不过因为不用乘原根，所以可以省去中间过程的取模，常数小得多。

子集卷积

注意到 i\ |\ j=k 就是 FWT_or，但是有一个 i\ \&\ j=0 的限制，相当于说 i,j 是 k 的一个划分。
考虑转化“与”条件，设 pc(i) 表示 i 二进制下 1 的数量，那么“与”条件可以变为 pc(i)+pc(j)=pc(k)。这非常像多项式乘法——我们给每一项 i 配一个 x^{pc(i)}，然后再把两个集合卷积即可。
（可以把原先数组的每一项都看作是一个多项式，理解起来可能会简单一点）

for(int i = 0; i <= n; i++) FWT_or(a[i], 1 << n, 1), FWT_or(b[i], 1 << n, 1);
for(int i = 0; i <= n; i++)
    for(int j = 0; j <= i; j++)
        for(int k = 0; k < 1 << n; k++)
            c[i][k] = (c[i][k] + 1ull * a[j][k] * b[i - j][k]) % mod;
for(int i = 0; i <= n; i++) FWT_or(c[i], 1 << n, -1);