题解：P6115 【模板】整式递推

· · 题解

首先考虑下面的问题：给定 d 次多项式 F(x)，求 \displaystyle\prod_{i=1}^nF(i)，这可以看作是本题 m=1 的特殊情况，例如取 F(x)=x 就是求 n!。

如果 F 是固定的，经典的方法就是分段打表，即设定块长 B，设 T=\left\lfloor\dfrac nB\right\rfloor，预先用 \mathcal O(nd) 的时间处理出 F!(B),F!(2B),\cdots,F!(TB) 并写到代码里，查询时再暴力处理散块即可，复杂度为 \mathcal O(Bd)，主要限制是代码长度，以及必须固定模数。

现在 F 是输入给定的，因此我们希望程序完成这个过程。设 G_k(x)=\displaystyle\prod_{i=1}^kF(x+i)，要求的就是 G_B(0),G_B(B),\cdots,G_B(TB)。直接倍增乘出 G_B(x) 然后多点求值的话是 \mathcal O\left(Bd\log n+\dfrac nB\log^2n\right) 的，取 B=\mathcal O\left(\sqrt{\dfrac{n\log n}{d}}\right) 有最优复杂度 \mathcal O(\sqrt{nd\log n}\log n)。

考虑一些更优秀的做法，比如说倍增，由于 \deg G_k=kd，我们至少需要 kd+1 个点值来表示 G_k，当然我们选择 G_k(0),G_k(B),\cdots,G_k(kdB)，这样最后刚好就能把答案求出来——这句话的前提是有 dB\ge T，即 dB^2\ge n。

为了将 k 倍增到 B，要做的就是将 k 加 1，以及将 k 乘 2。而加 1 是容易的，对于 i\le kd，有 G_{k+1}(iB)=G_k(iB)\cdot F(iB+k+1)，可以单个 \mathcal O(d) 更新；而对于 kd<i\le(k+1)d，直接单个 \mathcal O(kd) 暴力计算即可，这部分的复杂度为 \mathcal O(kd^2)。

考虑如何乘 2，有 G_{2k}(iB)=G_k(iB)\cdot G_k(iB+k)，因此需要额外得知 G_k((kd+1)B),G_k((kd+2)B),\cdots,G_k(2kdB) 和 G_k(k),G_k(B+k),\cdots,G_k(2kdB+k) 这些点值，前者相当于将原来的前 kd 个点值整体平移 (kd+1)B，后者相当于将得到的 2kd+1 个点值再整体平移 k。

现在要做的就是说，给定次数 \le n 的多项式 F(x) 的点值 F(0),F(l),\cdots,F(nl)，求 F(k),F(l+k),\cdots,F(nl+k)。根据拉格朗日插值公式，有：

\begin{aligned} F(il+k)&=\sum_{j=0}^nF(jl)\prod_{t\neq j}\dfrac{(i-t)l+k}{(j-t)l}=\sum_{j=0}^nF(jl)\prod_{t\neq j}\dfrac{i-t+k/l}{j-t}\\ &=(i+k/l)^{\underline{n+1}}\sum_{j=0}^n\dfrac{F(jl)}{(-1)^{n-j}j!(n-j)!}\cdot\dfrac{1}{i-j+k/l} \end{aligned}

前面的下降幂可以滑动窗口求，而后面的 \sum 是卷积的形式，只不过可能搞到负的下标去所以要平移一下；同时可能遇到一些乘除 0 的情况，仅当 l\mid k 时可能发生，且此时 F(il+k) 一定在一开始给定的 n+1 个点值之中，特判掉即可（更好的做法是，在传进去的时候就确保没有相交的点值）。因此将 k 乘 2 可以做到 \mathcal O(kd\log n)。

总复杂度为 T(k)=T(k/2)+\mathcal O(kd(d+\log n))，即 \mathcal O(Bd(d+\log n))，由于 dB^2\ge n，取 B=\sqrt{n/d} 有最优复杂度 \mathcal O(\sqrt{nd}(d+\log n))。

考虑本题：对于 n\ge m 有 \displaystyle\sum_{i=0}^m P_i(n)a_{n-i}=0，而 m,d 均较小，且 P_0(m)\sim P_0(n) 均不为 0，要求的是 a_n。

现在有 m>1 了，根据递推的通法，考虑构建转移矩阵：

\begin{bmatrix} -P_1(n)&-P_2(n)&-P_3(n)&\cdots&-P_m(n)\\ P_0(n)&0&0&\cdots&0\\ 0&P_0(n)&0&\cdots&0\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&0&\cdots&0 \end{bmatrix} \begin{bmatrix} a_{n-1}\\a_{n-2}\\a_{n-3}\\ \vdots \\a_{n-m} \end{bmatrix} = P_0(n) \begin{bmatrix} a_n\\a_{n-1}\\a_{n-2}\\ \vdots \\a_{n-m+1} \end{bmatrix}

有转移矩阵 M(n)=\left[\begin{array}{c|c}-\bm{p}(n)&-P_m(n)\\ \hline P_0(n) I&\bm{0}\end{array}\right]，只要求出 M(n)M(n-1)\cdots M(m+1)M(m) 即可得到 a_n。

依然考虑和上面类似的方法，设 L_k(x)=M(x+k)\cdots M(x+2)M(x+1) 也是 m\times m 的矩阵，里面每个元素都是 kd 次的多项式，点值 L_k(0),L_k(B),\cdots,L_k(kdB) 也是矩阵，而乘法是矩阵乘法，考虑如何修改之前的做法。

其实大差不差，将 k 加 1 的部分直接照抄即可，复杂度是 \mathcal O(kd^2(m^2d+m^3)) 的，瓶颈在暴力算 (kd,(k+1)d] 内的 i，应该可以改成点值平移，不过大概率会更慢；将 k 乘 2 的部分，问题在于点值平移时会有矩阵序列卷上整数序列的情况，但显然平移时矩阵中每个元素都是互相独立的，把 m^2 个位置单独拿出来做即可，复杂度为 \mathcal O(kd(m^2\log n+m^3))；取 B=\sqrt{n/d}，有最优复杂度 \mathcal O(\sqrt{nd}(m^2d^2+m^3d+m^2\log n))。

这也太抽象了，不过可以偷个懒，把 B 取成 \ge\sqrt{n/d} 的最小的 2 的整数次幂，也就是说只用写将 k 乘 2 的部分，这样复杂度就是 \mathcal O(\sqrt{nd}(m^3+m^2\log n)) 的了，不过这样比较需要注意常数，一种卡常的方式是，多项式点值平移的时候，是一对下标范围 [0,n] 和 [0,2n] 的数组做卷积，正常做需要 3n+1 的卷积长度；不过最后只需要 [n,2n] 内的值，由于 NTT 是循环卷积，取 [0,2n) 的卷积范围，这样下标 n 位置会被原本的 3n 位置影响，且 2n 位置加到了 0 位置上，额外计算贡献即可；由于 k 是 2 的整数次幂，这样就能少一半的 NTT 长度，但实现会变得很复杂。

实现时需要特别注意矩阵乘法的左右顺序，例如倍增时有 L_{2k}(iB)=L_k(iB+k)\cdot L_k(iB)；别忘了最后还要除掉 \displaystyle\prod_{i=m}^nP_0(i)，直接按照一开始的办法求出即可。

还能不能做到更优呢？感觉是很困难的。

代码写的不太好看就不放了，而且这个题 AC 记录怎么有一大半都是超的题解，害怕。