整式递推机械化

qwaszx · 2021-01-25 09:01:09 · 个人记录

写了两天的成果，作为娱乐以及工业代码的练习

基于 EI 的 ODE 自动机

加法和乘法的处理和 EI 是一样的，不过复合不是很一样。

现在我们知道 f 满足的 ODE f^{(n)}=\sum_ip_if^{(i)} 以及代数方程 P(u)=0，其系数都 \in K(x)，即为 K 上的有理分式，我们想要求出一个 f\circ u 满足的 ODE.

为了简化问题，只考虑 P(u) 为不可约多项式（指不能被分解为两个更低次的系数 \in K(x) 的多项式的积）的情况，这也和我们的实际应用相符.

先考虑简单的情况，即 u 为有理分式的情况（或者说 P(u) 为一次的情况），这时我们考虑

(f\circ u)^{(k)}=\sum_{i}a_{k,i}f^{(i)}\circ u

对其求导时，我们考虑某一项的导数

(a_{k,i}(f^{(i)}\circ u))'=a_{k,i}'(f^{(i)}\circ u)+a_{k,i}(f^{(i+1)}\circ u)u'

当 i 足够大时我们就可以利用 f 的 ODE 来降次，降次时需要做一次系数的复合:

f^{(n+1)}\circ u=\sum_i(p_i\circ u)(f^{(i)}\circ u)

接下来考虑 P(u) 高于一次的情况，这时我们需要处理两个问题.

把 u' 用 u 表示

我们有
\frac{\partial P(u,x)}{\partial x}=\frac{\partial P}{\partial u}u'+\frac{\partial P}{\partial x}=0
于是就有
u'=\frac{-\partial_xP}{\partial_uP}
接下来要把 u' 化为关于 u 的整式，我们可以将 u' 对 P(u) 取模，这就要求一个同余方程 k\partial_uP\equiv -\partial_x P\pmod{P(u)}. 这是否一定有解呢？答案是肯定的，我们只需要证明 \gcd(P(u),\partial_u P)|\partial_x P，为此我们考虑 P(u) 的因式分解：
P(u)=c\prod_i p_i(u)^{k_i}
那么我们知道
\gcd(P(u),\partial_u P)=\prod_i(p_i(u))^{k_i-1}
而
\partial_x P=\sum_i k_ip_i(u)^{k_i-1}(\partial_x p_i)\prod_{j\neq i}(u-a_j)^{k_j}
于是显然任何一个 p_i(u)^{k_i-1} 都整除 \partial_x P.

于是我们只需要用多项式扩展欧几里得算法（O(n^2) 即可）计算出这个方程的解即可把 u' 表示为 u 的整式.
计算 p_k\circ u

问题还是一样的，我们要把一个分式对 P(u) 取模从而化成整式. 由于我们的需要，这个分式的系数都是 K 上的数，而 P(u) 的根应该都是在 K(x)\backslash K 上的，于是它们没有公因式，逆元存在.

到这里原理就阐述清楚了，不过实现上还是有些挺有意思的地方

我们首先定义一个模板类 Poly<T>，表示一个 T 上的多项式. 由于长度总是很小，运算使用 n^2 暴力计算.
接下来定义模板类 FRAC<T>，表示一个 T 上的有理分式.
为了方便，还要定义模域 Z，表示域 Z_{\mathrm{mod}}，用于实现通常的取模意义下的操作.
最后我们定义模板类 ODE<T>，表示一个系数在 T(x) 上的 ODE.

那么通常的多项式板子可以认为是特化了 Poly<Z> 的情况. 这套体系有什么方便之处呢？答案是我们可以更方便地进行符号计算，例如我们经常需要计算 x^k\circ u 的 ODE，而 k 不是固定的. 那么我们可以使用 ODE<FRAC<Z>> 来完成系数带着 k 的计算. 因为 ODE 自动机直接计算比较慢（当然还是比 n\log n 快得多），所以如果要达成更高的效率需要得到带着 k 的整式递推之后特化实现.

这里是一个实现~~当然我觉得也没有人会看的~~