高斯消元没 有 精 神，我们暴力维护分式基才是你们的老 大 哥！
看到有些毒瘤把中国象棋这道题整出来一个彪悍的递推式但又很可惜代码不是线性，于是就写了一个让电脑帮我自动算整式递推的非 BM 类方法。
我们先回顾一下中国象棋这道题，但是题目本身不是本篇文章的重点，我们先快进到生成函数
∑ n a n s n , n + k x n n ! ( n + k ) ! = ( 1 − x ) − 1 / 2 exp ( x ( 1 + x ) 2 − 2 x ) ( 2 − x 2 − 2 x ) k k ! ⋅ 0 F 1 ( ; k + 1 ; x ( 2 − x 2 − 2 x ) 2 ) \sum_n \mathrm{ans}_{n,n+k} \frac{x^n}{n!(n+k)!} = (1-x)^{-1/2} \exp\left(\frac{x(1+x)}{2-2x}\right) \frac {\left(\dfrac{2-x}{2-2x}\right)^k}{k!} \cdot { {}_0 F_1\left(;k+1; x\left(\frac{2-x}{2-2x}\right)^2\right) } n∑ansn,n+kn!(n+k)!xn=(1−x)−1/2exp(2−2xx(1+x))k!(2−2x2−x)k⋅0F1(;k+1;x(2−2x2−x)2)
为了最大幅度的简化我们人类的脑力劳动，我们首先来看看我们可以尽量使用哪些法则来简化我们需要做的工作。
我们首先确认几个简单的事实：
设 f = x α f=x^\alpha f=xα，可得 f ′ = α x α − 1 f'=\alpha x^{\alpha -1} f′=αxα−1，既得 α f − x f ′ = 0 \alpha f-xf'=0 αf−xf′=0。
指数函数 f = exp x f=\exp x f=expx 满足 f − f ′ = 0 f-f'=0 f−f′=0，不必多说。
定义算子 ϑ = x D \vartheta = xD ϑ=xD，广义超几何函数
f = p F q ( a 1 , … a p ; b 1 , … , b q ; x ) = ∑ n a 1 n ‾ ⋯ a p n ‾ b 1 n ‾ ⋯ b q n ‾ n ! x n f = {}_p F_q(a_1,\dots a_p; b_1, \dots, b_q;x) = \sum_n \frac{a_1^{\overline n} \cdots a_p^{\overline n}}{b_1^{\overline n} \cdots b_q^{\overline n}n!}x^n f=pFq(a1,…ap;b1,…,bq;x)=n∑b1n⋯bqnn!a1n⋯apnxn
满足
( ϑ + a 1 ) ⋯ ( ϑ + a p ) f = D ( ϑ + b 1 − 1 ) ⋯ ( ϑ + b q − 1 ) f (\vartheta + a_1) \cdots (\vartheta + a_p)f = D (\vartheta + b_1 - 1) \cdots (\vartheta + b_q - 1)f (ϑ+a1)⋯(ϑ+ap)f=D(ϑ+b1−1)⋯(ϑ+bq−1)f
把两边的算子展开，就都是 ODE 了。
并且我们能够从书中查阅到一些如下结论：
若 f , g f,g f,g 微分有限，那么 f + g , f g f+g,fg f+g,fg 微分有限。若 f f f 微分有限， g g g 是代数的，那么 f ∘ g f\circ g f∘g 微分有限。那么我们发现基于以上法则，我们已经完全能判断出答案是整式递推了。
其中的证明基本上是围绕 f , f ′ , … f, f', \dots f,f′,… 构成的线性空间考虑的，这其实就启发了我们如何通过 f , g f,g f,g 的 ODE 维护出他们组合起来所满足的 ODE。
给定 f , g f,g f,g 分别满足的 ODE，我们希望找出一个 f + g f+g f+g 必然满足的 ODE。我们不妨设
f ( n ) = ∑ 0 ≤ k < n p k f ( k ) f^{(n)} = \sum_{0\le k<n} p_k f^{(k)} f(n)=0≤k<n∑pkf(k)
其中 p k ∈ K ( x ) p_k \in K(x) pk∈K(x) 即是 K K K 上的有理分式。类似的， g ( m ) g^{(m)} g(m) 有类似的系数 q k q_k qk。
接下来我们考虑如下函数列 ( f + g ) , ( f + g ) ′ , ( f + g ) ( 2 ) , … (f + g), (f+g)', (f+g)^{(2)},\dots (f+g),(f+g)′,(f+g)(2),…：
不妨设 n ≤ m n\le m n≤m，当 k < n k<n k<n 的时候， ( f + g ) ( k ) = f ( k ) + g ( k ) (f+g)^{(k)} = f^{(k)} + g^{(k)} (f+g)(k)=f(k)+g(k) 。而当 k = n k=n k=n 的时候，我们注意到 f ( n ) f^{(n)} f(n) 就可以通过我们满足的微分方程来降次了。因此我们总是可以将 ( f + g ) ( k ) (f+g)^{(k)} (f+g)(k) 表示为
( f + g ) ( k ) = ∑ i < n a i f ( i ) + ∑ j b j g ( j ) (f+g)^{(k)} = \sum_{i<n} a_i f^{(i)} + \sum_j b_jg^{(j)} (f+g)(k)=i<n∑aif(i)+j∑bjg(j)
进而有
( f + g ) ( k + 1 ) = ( ∑ i < n a i f ( i ) + ∑ j b j g ( j ) ) ′ = ∑ i < n ( a i ′ f ( i ) + a i f ( i + 1 ) ) + ∑ j ( b j ′ g ( j ) + b j g ( j + 1 ) ) \begin{aligned} (f+g)^{(k+1)} &= \left(\sum_{i<n} a_i f^{(i)} + \sum_j b_jg^{(j)}\right)'\\ &= \sum_{i<n} \left(a_i' f^{(i)} + a_i f^{(i+1)} \right) + \sum_j \left( b_j'g^{(j)} + b_j g^{(j+1)}\right) \end{aligned} (f+g)(k+1)=(i<n∑aif(i)+j∑bjg(j))′=i<n∑(ai′f(i)+aif(i+1))+j∑(bj′g(j)+bjg(j+1))
因此我们考虑向量 ( a 0 ⋯ a n − 1 b 0 ⋯ b m − 1 ) T \begin{pmatrix}a_0 & \cdots & a_{n-1} & b_0 & \cdots & b_{m-1}\end{pmatrix}^{\mathsf T} (a0⋯an−1b0⋯bm−1)T 本来就只有 n + m n+m n+m 个维度，因此 ( u + v ) ( 0 ) , … , ( u + v ) ( n + m ) (u+v)^{(0)}, \dots, (u+v)^{(n+m)} (u+v)(0),…,(u+v)(n+m) 必然线性相关。我们只需要像维护线性基一样不断的把 ( u + v ) ( k ) (u+v)^{(k)} (u+v)(k) 其被表示的向量丢进去，直到有一个向量被以前的向量线性表出即可。这也证明了 f + g f+g f+g 满足的 ODE 的微分次数是不超过 n + m n+m n+m 的。
我们类似地考虑
( f g ) ( k ) = ∑ i < n , j < m a i , j f ( i ) g ( j ) (fg)^{(k)} = \sum_{i<n,j<m} a_{i,j} f^{(i)}g^{(j)} (fg)(k)=i<n,j<m∑ai,jf(i)g(j)
这种表示方法，它在微分的时候满足 ( t f ( i ) g ( j ) ) ′ = t ′ f ( i ) g ( j ) + t f ( i + 1 ) g ( j ) + t f ( i ) g ( j + 1 ) (tf^{(i)}g^{(j)})'=t'f^{(i)}g^{(j)} + tf^{(i+1)}g^{(j)} + tf^{(i)}g^{(j+1)} (tf(i)g(j))′=t′f(i)g(j)+tf(i+1)g(j)+tf(i)g(j+1)。
因此只有 n m nm nm 个维度，我们类似地维护即可。
对于一般性的 g ∈ K a l g [ [ x ] ] g\in K_{\mathrm{alg}}[[x]] g∈Kalg[[x]] 即在 K [ x ] K[x] K[x] 上代数的情况比较麻烦，这里给出 g ∈ K ( x ) g\in K(x) g∈K(x) 的情况，即有理分式。
我们考虑链式法则 ( f ∘ g ) ′ = ( f ′ ∘ g ) ⋅ g ′ (f\circ g)' = (f'\circ g) \cdot g' (f∘g)′=(f′∘g)⋅g′，根据分式求导发展显然 g ′ ∈ K ( x ) g' \in K(x) g′∈K(x)，我们可以得到一个下三角矩阵
( f ∘ g ) ( k ) = ∑ j a k , j ( f ( k ) ∘ g ) (f\circ g)^{(k)} = \sum_j a_{k,j} (f^{(k)} \circ g) (f∘g)(k)=j∑ak,j(f(k)∘g)
设向量 v \mathbf v v 的分量 v k v_k vk 是 ODE 中 f f f 的 k k k 阶导上的系数。不难发现我们所得到的新 ODE 就是
( v ∘ g ) A − 1 h = 0 (\mathbf v \circ g) \mathbf A^{-1} \mathbf h = 0 (v∘g)A−1h=0
其中 h k = ( f ∘ g ) ( k ) h_k = (f\circ g)^{(k)} hk=(f∘g)(k)。
我们实际上需要计算的就是 ( v ∘ g ) A − 1 (\mathbf v \circ g) \mathbf A^{-1} (v∘g)A−1，这个只需对增广矩阵 [ ( v ∘ g ) ∣ A ] [(\mathbf v \circ g) | \mathbf A] [(v∘g)∣A] 做消元，把右边消成 I \mathbf I I 的时候就得到了。注意因为是上三角矩阵，所以只需要 Θ ( n 2 ) \Theta(n^2) Θ(n2) 次实际上是单点的行列变换。
这个东西看起来很有意思，但是实际上有一个诡异的点，就是在求 f + g , f g f+g, fg f+g,fg 这种的过程中会发生严重的 中间表示膨胀 (intermediate expression swell) 问题，因为我们在反复求导的过程中系数有分式，而分式求导中 ( u v ) ′ = u ′ v − v ′ u v 2 \left(\frac uv\right)'=\frac{u'v-v'u}{v^2} (vu)′=v2u′v−v′u 我们发现分母的次数是倍增的，因此在过程中分子可能会以指数级增长。经过测试，代码在寻找 ( ln ( 1 − x ) ) 3 (\ln(1-x))^3 (ln(1−x))3 的 ODE 的时候在本机花了 3ms，中间通分过程最高产生过一个 84 次多项式，尽管结果的多项式并不大。
这里是一份 chess AC code，大家也可以玩一玩试试看，并且一些注释没有删，可以通过信息看到这个是怎么执行的。
