数学杂谈 #21
初探整式递推和 ODE
约定
记 \(D=\frac{\text d}{\text dx},\vartheta=xD\)。前者是微分算子,后者可以看作对于形式幂级数,将每一项的指数乘到它的系数里去。
记 \(f^{(k)}=D^kf\)。
ODE
ODE 和 D-finite
ODE 的意思是“常微分方程(Ordinary Differential Equation)”。
Note.
查了一下百科,Ordinary 指的是“只有一个变元”,也就是未知函数是一元函数的微分方程。
它的一般形式是:
其中 \(n\) 是某个正整数。(据说)OI 中常用的形式为:
(据说)\(p_k(x)\) 一般是常数次数或者常数项数的有理分式,\(n\) 一般是常数。特别地,我们要求 \(p_n(x)\neq 0\),否则这个微分方程不是阶数最小的。
对于一个函数 \(f(x)\),若它满足一个上述常微分方程,那么就称它为“微分有限(D-finite)”的。
常见的 D-finite 的函数
有限多项式
如果 \(A(x)=\sum_{k=0}^ma_kx^k\),显然有 \(A^{(m+1)}=0\)。
这个东西没有 \(\Delta^{m+1}A=0\) 好用 qwq。
当然,如果差分对于下降幂来说是平凡的位移效果,而在普通幂上面可以建立方程,我们也许可以考虑微分作用在下降幂上的结果。
幂函数
如果 \(f(x)=x^a,a\neq 0\),显然有 \(af-xf'=0\)。
指数函数
如果 \(f(x)=e^x\),显然有 \(f-f'=0\)。其它的 \(f=a^x\) 可以通过复合来表达。
对数函数
如果 \(f(x)=\ln x\),显然有 \(xf'-1=0\)。其它的 \(f=\log_ax\) 可以通过复合来表达。
注:虽然这个东西不是一般的线性组合形式,但是它仍然满足一般的 ODE 的形式。
三角函数和反三角函数
三角的 ODE 很简单,但是没啥用,不写了。
反三角的 ODE 可以通过查导数得到,并且更加没用,也不写了。
超几何级数
注:以下内容均来自于《具体数学》和
slide.pdf
,我都是抄的。
超几何级数的定义是:
其中任意的 \(b\) 都不可以是非正整数。
根据 slide.pdf
,我们可以得到超几何级数的一般的 ODE 为:
证明不会。但是超几何级数相当之强悍,它涵盖了超级多的常见的幂级数。比如几何级数 \(\frac{1}{1-x}\)、指数函数 \(e^x\)、(几乎是)对数函数 \(x^{-1}\ln(1+x)\) 等等。所以记一下也无妨。
Note.
多说几句。《具体数学》上说,超几何级数涵盖了所有“相邻两项的系数之比是 \(k\) 的有理分式”的级数。简单来说,如果令 \(c_k\) 为:
\[\begin{aligned} c_k=\frac{\prod_{j=1}^ma_{j}^{\overline k}}{k!\prod_{t=1}^nb_{t}^{\overline k}} \end{aligned} \]比一下就是:
\[\frac{c_{k+1}}{c_k}=\frac{k!}{(k+1)!}\frac{\prod_{j=1}^ma_j^{\overline{k+1}}}{\prod_{j=1}^ma_j^{\overline{k}}}\frac{\prod_{t=1}^nb_t^{\overline{k}}}{\prod_{t=1}^nb_t^{\overline{k+1}}}=\frac{\prod_{j=1}^m(a_j+k)}{(k+1)\prod_{t=1}^n(b_t+k)} \]
ODE 与函数运算
我们已经得到了五类基本函数的 ODE,之后只需要理解在若干函数运算下 ODE 的相应变化,我们(理论上)就可以给出任意初等函数的 ODE。
在此之前,我们先来简单看一条比较重要的引理:
Lemma.
形式幂级数 \(u(x)\) 微分有限当且仅当 \(\dim_{\mathbb K(x)}\operatorname{span}_{\mathbb K(x)}\{u,u',u'',\dots\}\) 有限。
从右到左是极其容易的。从左到右感性理解就好,如果要深入理解证明思路,看一看下面“加法”部分就明白了。
加法
Theorem.
如果 \(f,g\) 是微分有限的,则 \(f+g\) 也是微分有限的。
Proof.
假设 \(f,g\) 各自满足微分方程:
\[\sum_{k=0}^np_k(x)f^{(k)}=0\\ \sum_{k=0}^mq_k(x)g^{(k)}=0 \]则我们只需要构建一个关于 \((f+g)\) 的微分方程即可完成证明。
方程本身给出了关系:
\[f^{(n)}=-p_n(x)^{-1}\sum_{k=0}^{n-1}p_k(x)f^{(k)} \]两边同时对 \(x\) 求导并用上式替换 \(f^{(n)}\),我们即可得到 \(f^{(n+1)}\) 关于 \(f,f',f'',\dots,f^{(n-1)}\) 的线性表示。
重复这个过程,我们可以得到 \(f^{(n+k)}\) 关于 \(f,f',f'',\dots,f^{(n-1)}\) 的线性表示,其中 \(k\in \mathbb N\)。
\(g\) 类似。于是任意的 \((f+g)^{(k)}\) 都可以被 \(f,f',f'',\dots,f^{(n-1)},g,g',g'',\dots,g^{(m-1)}\) 线性表示。
从而 \((f+g),(f+g)',(f+g)^{''},\dots,(f+g)^{(n+m)}\) 必然线性相关,所以存在一个阶数不超过 \(n+m\) 的 ODE。实现的时候只需要维护 \((f+g)^{(k)}\) 的线性基即可。
乘法
Theorem.
如果 \(f,g\) 是微分有限的,则 \(fg\) 也是微分有限的。
Proof.
使用公式:
\[(fg)^{(n)}=\sum_{k=0}^{n}\binom{n}{k}f^{(k)}g^{(n-k)} \]所以任意的 \((fg)^{(n)}\) 均可被 \(f^{(k)}g^{(j)}\) 线性表示,其中 \(0\le k<n,0\le j<m\)。故必然存在一个阶数不超过 \(nm\) 的 \(fg\) 的 ODE。
复合
Definition.
一个形式幂级数 \(A(x)\) 是“代数形式幂级数”,当且仅当 \(A(x)\) 可以作为某个系数为 \(\mathbb K(x)\) 的多项式方程的根。其中 \(\mathbb K\) 为某个域,\(\mathbb K(x)\) 表示系数在 \(\mathbb K\) 内的有理分式域。
我们将要指出:
Theorem.
如果 \(f\) 微分有限,而 \(g\) 是一个代数形式幂级数,则 \(f\circ g\) 也微分有限。
Explanation.
考虑 \(g\) 是有理分式的情况。反复使用链式法则我们可以得到:
\[(f\circ g)^{(k)}=\sum_{j=0}^{k}c_{k,j}(x)\left(f^{(j)}\circ g\right) \]\(c_j(x)\) 是关于 \(g\) 及其多阶微分的线性组合,但是 \(g\) 是有理分式。
而 \(f^{(j)}\circ g\) 的 ODE 可以直接由 \(f\) 的 ODE 生成:
\[\sum_{k=0}^np_k(g)(f^{(k)}\circ g)=0 \]因此 \(f^{(k)}\circ g\) 构成一组基,其中 \(0\le k<n\)。进而必然存在一个阶数不超过 \(n\) 的 ODE。
假设 \(g\) 满足代数方程 \(P(x,g)=0\),其中 \(P(x,g)=\sum_{k=0}^{m-1}b_k{(x)}g^k\)。
仍然试图得到使用链式法则,但是这个时候我们就必须直面 \(g\) 的导数了。
先考虑 \(g'\),我们来通过 \(P(x,g)\) 生成关系。在方程两侧同时对于 \(x\) 求导:
\[\frac{\text d}{\text dx}P(x,g)=\frac{\partial}{\partial x}P(x,g)+\frac{\partial}{\partial g}P(x,g)\cdot g'=0 \]所以有 \(g'=-\frac{\frac{\partial}{\partial x}P(x,g)}{\frac{\partial}{\partial g}P(x,g)}\)。上下都是关于 \(g\) 的多项式,我们尝试将 \(g'\) 也写成关于 \(g\) 的多项式。
我们知道 \(P(x,g)=0\),因此对于 \(P(x,g)\) 做截取。也即,我们只需要求:
\[\frac{\partial}{\partial x}P(x,g)+\frac{\partial}{\partial g}P(x,g)u(x,g)\equiv 0\pmod{P(x,g)} \]“可以证明这个总有解,为了防止多解,我们再假设 \(P(x,g)\) 关于 \(g\) 在 \(\mathbb K(x)\) 上不可约”。解出来 \(u(x,g)\) 即可。据说这个过程需要多项式欧几里德,不过考虑到平时大家都用手算,无所谓了。😄
另一方面,我们回头看一下:
\[\sum_{k=0}^np_k(g)(f^{(k)}\circ g)=0 \]怎么算 \(a_k(g)\)?“可以使用多项式欧几里德”。我觉得大家手算的时候还是直接复合好了。
如果 \(a_k(g)\) 的分子分母都比较短,我们把所有 \(a_k\) 通分,之后去掉分母即可将有理分式变成有限多项式。
最后补充说明一下(不是泼冷水),下面很多 ODE 都可以通过反复求导建立,所以不要忘了这个基本策略!
整式递推
整式递推是一类并不特殊的递推数列:
Definition.
对于无穷数列 \(\{a_k\}_{k\ge 0}\) 和有限非空多项式列 \(\{r_0,r_1,r_2,\dots,r_{m}\}\),如果对于任意的 \(n\ge m\),都有:
\[\sum_{k=0}^ma_{n-k}r_{k}(n)=0 \]则我们称 \(r\) 为 \(a\) 的整式递推式。我们称存在整式递推式的无限数列为整式递推数列。
对于有穷数列 \(\{a_k\}_{k=0}^{p}\),我们可以类似地定义。此处略去不谈。
整式递推涵盖的范围很广。比如,如果所有的 \(r\) 均取到常数,则这个递推式退化成了线性递推。
下面这个定理很好地指出了整式递推和数列的 GF \(A(x)\) 的 ODE 的关系(之一):
Theorem.
一个数列 \(\{a_k\}_{k\ge 0}\) 是整式递推数列,当且仅当数列的 GF \(A(x)\) 是 D-finite 的。
显然,我们只需要给出一个构造性证明即可在证明定理的同时,给出从“整式递推”到“ODE”的转化途径。
从 ODE 到整式递推是极其容易的——我们只需要对于 ODE 取它的 \([x^n]\) 即可建立一个方程,怎么怎么移位一下就可以得到递推式了。
从整式递推到 ODE 反而不太容易。考虑到这一部分内容在 OI 中的应用基本都是前者,所以我并不准备展开讲。
例题
短分式幂
给定 \(F(x)=\frac{P(x)}{Q(x)}\),满足 \(P(x),Q(x)\) 均为多项式且 \(\max\{\deg P,\deg Q\}=m\)。现在请求出 \(\ln F(x),F(x)^\alpha,\exp F(x)\) 的前 \(n\) 项。你可以认为 \(m\ll n\),所以这里要求你给出一个 \(O(nm)\) 的算法。运算在某个域上进行。
经典问题了。
首先解决 \(\ln F(x)\)。显然有 \(\ln F(x)=\ln P(x)-\ln Q(x)\),所以我们这里仅考察 \(U(x)=\ln P(x)\) 的情况。
记得上面我们提到的 \(\ln x\) 的微分方程吗?这里我们还暂时用不到它(笑。不过我们可以以此获得启发,在等式两侧同时微分:
这样我们实际上就得到了一个 \(U(x)\) 的微分方程。作为例子,我们再详细一点——两侧同时取 \([x^n]\):
整理一下即可得到:
这相当于是建立了 \(\{u\}\) 的一个整式递推数列。由于 \(P(x)\) 比较短,因此 \(k\) 只需要枚举 \(O(m)\) 项。递推的复杂度就是 \(O(nm)\) 的。
其它两类照猫画虎。对于 \(V(x)=F(x)^\alpha\),两侧同时微分:
这里 \(\max\{\deg P'(x)Q(x),\deg Q'(x)P(x), \deg P(x)Q(x)\}\le 2m\)。因此我们可以先 \(O(m^2)\) 暴力卷积,之后再做递推。
对于 \(W(x)=\exp F(x)\),两侧同时微分:
类似处理即可。
我们在上面伏了一笔。上面的微分方程都是由普通的微分建立的,如何用上之前 ODE 与运算的结论呢?
举个例子,设 \(g(x)=x^\alpha\),那么 \(V=g\circ F\)。我们知道 \(g(x)\) 满足 \(\alpha g-xg'=0\),所以有 \(g'=\frac{\alpha g}{x}\)。
怎么变成 \(V\) 的微分方程?根据结论,\(V\) 及其任意阶微分必然都可以由 \(V=g\circ F\) 自身来线性表示。所以我们单独考察 \(V'\) 即可生成 ODE。
链式法则告诉我们 \(V'=(g\circ F)'=F'\cdot (g'\circ F)=F'\cdot \frac{\alpha V}{F}\),因此有 \(F(x)V'(x)=\alpha F'(x)V(x)\)。这和我们直接微分的结果(当然)是一致的。
「UOJ514」通用测评号
开门见山,我们需要求:
其中 \(A(x)=\sum_{k=0}^{a-1}\frac{x^kn^{-k}}{k!},B(x)=\sum_{k=0}^{b-1}\frac{x^kn^{-k}}{k!}\)。
数据满足 \(1\le n\le 250,1\le b<a\le 250\),你只需要给出一个 \(O(n^3)\) 的算法(认为 \(n,a,b\) 同阶)。
首先发现把式子内层不完全展开之后可以得到 \(x^{\alpha}e^{\beta xn^{-1}}\) 的线性组合。而 \(\sum_{k\ge 0}k![x^k]x^{\alpha}e^{\beta xn^{-1}}\) 是容易计算的。因此,我们的目标就是算出来每一个单项之前的系数。
使用二项式定理展开 \((e^{xn^{-1}}-B(x))^{n-2}\),我们相当于只需要快速算出 \(U_k(x)=A(x)B(x)^k,V_k(x)=B(x)^k\) 的所有系数。考虑到它们的项数之和就是 \(O(n^3)\) 级别的,我们能想到的一种比较好的方法,就是递推系数。
以较复杂的 \(U(x)\) 为例,我们在两侧同时微分:
注意到 \(A'(x),B'(x)\) 都是对于 \(e^{xn^{-1}}\) 进行截断之后的结果,这直接导致它们的 ODE 极其简单:
之后我们会详细讨论一下截断对于 ODE 的影响。
下面我们就记 \(\iota_A(x)=A(x)-A'(x),\iota_B(x)=B(x)-B'(x)\)。所以,我们可以将 ODE 改写成:
此外,我们显然还有边界值 \(U_0=A,V_0=B\)。因此,我们可以按照 \(k\) 从小到大递推,复杂度即为 \(O(n^3)\)。
「CTS2019」珍珠
有 \(n\) 个在范围 \([1,D]\) 内的整数均匀随机变量。
求至少能选出 \(m\) 个瓶子,使得存在一种方案,选择一些变量,并把选出来的每一个变量放到一个瓶子中,满足每个瓶子都恰好装两个值相同的变量的概率。
请输出概率乘上 \(D^n\) 后对 \(998244353\) 取模的值。
对于 \(100\%\) 的数据,满足 \(1\le D\le 10^5,0\le m\le 10^9,1\le n\le 10^9\)。
注意到,我们实际上只需要讨论一下有多少种值出现了奇数次即可。假如有 \(k\) 种值出现了奇数次,则我们要求 \(\lfloor\frac{n-k}{2}\rfloor\ge m\) 和 \(k\equiv n\pmod 2\)。
当我们写出 EGF 的时候,后面的同余限制实际上不用管,因为不满足的项根本不会贡献。假设第一条给出限制为 \(k\le K\),答案为:
注意到和式内的其实是 \(e^x\) 的生成函数。因此我们可以换元得到 \(F(z)\):
现在我们就可以直接计算 \(F(z)\) 了,据说可以 \(O(D\log^2 D)\) 算。如果将和式展开我们还可以得到 \(O(D\log D)\) 的做法。
知道了 \(F(z)=\sum_j f_jz^j\) 之后,我们相当于要计算 \(n!2^{-D}\sum_j f_j[x^n]e^{jx}\),这个就可以直接计算了。
但是这样还不够好。我们将 \(F(z)\) 处理成:
再次换元。设 \(H(z)=\frac{z^2-1}{z^2+1},G(t)=\sum_{k=0}^K\binom{D}{k}t^k\),现在就有 \(F(z)=z^{-D}(z^2+1)^DG(H(z))\)。
压力来到了 \((z^2+1)^DG(H(z))\) 身上,我们尝试给出这样一个 GF 的 ODE,这样就可以递推了。
先从最简单的 \(G(t)\) 开始。我们发现有 \(G(t)=(1+t)^D\bmod t^{K+1}\)。这能带来什么信息呢?
没有截断的 ODE 我们是知道的。那么,引入截断的 ODE 必然可以表示成:
这个 \(D(t)\) 的性状取决于 LHS 的抵消情况。考察 LHS 的某一项系数:如果它的系数全由 \(G\) 中(在微分前的)原先指数 \(\le K\) 的项提供,那么这一项系数不会受到截断的影响;如果它的系数全由 \(G\) 中(在微分前的)原先指数 \(>K\) 的项提供,那么这一项系数必然全是 \(0\)。所以 \(D(t)\) 包含的就是 LHS 中系数同时由 \(\le K\) 和 \(>K\) 的项提供的那些项。
依据这样的分析,我们可以发现 \(\mathscr D(t)=(K+1)\binom{D}{K+1}t^K\),因为 \(1\times G'(t)\) 恰好在 \(t^K\) 的位置缺了一项。因此:
现在我们可以快速地递推 \(G\) 的系数了。
下一步,考虑 \(G(H(z))\) 的 ODE。从这里开始,我们设 \(P(z)=G(H(z))\)。
然而我不想动脑筋了,直接使用 ODE 机械化构造的方法。
首先,给出 \(G'(t)=\frac{DG(t)-\mathscr D(t)}{1+t}\)。根据结论,我们考察 \((G\circ H)'\) 即可得到一个 ODE。有 \((G\circ H)'=H'\cdot (G'\circ H)\),所以有:
这个其实也可以线性推系数。
手算得到 \(H'(z)=\frac{4z}{(z^2+1)^2},1+H(z)=\frac{2z^2}{z^2+1}\),二者恰好抵消得到 \(\frac{1}{2}z(z+1)\),因此两个系数都很短。
\(H(z)^K\) 的系数虽然不方便直接算出来,但是我们可以使用短分式幂的方法和 \(P(z)\) 同步递推。
最后,我们只需要给出 \((z^2+1)^DP(z)\) 的 ODE 就好啦!仍然使用结论,\((z^2+1)^D\) 和 \(P(z)\) 各自都有一个一阶 ODE,因此表示 \(F(z)\) 及其任意次微分的结果也只需要用 \((z^2+1)^D\cdot P(z)\) 一项。因此直接考察 \(F'(z)\) 肯定可以得到一个 ODE:
整理一下即可得到最终的 ODE:
递推系数都比较短。后面的 \(R(z)=(z^2-1)^K(z^2+1)^{D+1-K}\) 可以建立另外的微分方程:
同步递推即可。这样我们就可以正儿八经地做到 \(O(D)\) 了。
初探载谭
先让我们从一道例题开始。
「CF932E」Binomial Sum
给定 \(n,m\),求出:
请给出一个 \(O(m+\log n)\) 的算法。
我们可以敏锐地觉察到,这个问题要求的其实就是:
这样通过基础多项式运算我们可以得到一个 \(O(m\log m)\) 的算法,已经相当不错了。
我们依然可以注意到,这个 GF 是一个 \(F(G(x))\) 的复合形式,其中 \(F(z)=(1+z)^n,G(x)=e^x\)。
回忆一下,复合有两种良好定义:
对于 \(F(G(x))\),当 \(F\) 有限或者 \([x^0]G(x)=0\) 的时候,该复合是良定义的。
这里虽然 \((1+x)^n\) 已经满足了“\(F\) 有限”的条件,但是它太长长长长长了,不方便我们处理。我们尝试将内层函数变成 \([x^0]G(x)=0\) 的形式。
这很简单,由于 \(G(0)=1\),我们只需要在 \(1\) 处展开 \(F(z)\) 即可:
关键的一步来了:由于 \([x^0](e^x-1)=0\),所以我们的求和上界只需要保留到 \(m\)。也就是说,当我们求 \([x^m]F(G(x))\) 的时候,我们用 \(\mathscr F(z)=F(z)\bmod (z-1)^{m+1}\) 不会导致结果出错,还可以降低次数!
现在,我们可以先算出来 \(F(z+1)\bmod z^{m+1}\),然后做一个多项式平移得到 \(\mathscr F(z)\)。这个算法仍然是 \(O(m\log m)\) 的。
但是 \(F(z),\mathscr F(z)\) 的形式都很简洁,我们尝试给出一个 ODE。先给出 \(\mathscr F(z+1)\) 的 ODE——我们之前已经研究过的那一个:
现在再平移就 OK 辣!
我们还需要解决 \([z^0]\mathscr F(z)\) 的边界情况。实际上只需要带进去算一算:
就好了。之后系数是可以 \(O(m)\) 递推的。而答案就是 \(m![x^m]\mathscr F(e^x)\),线性筛整数幂就好了。
EI 的载谭
我们把上面的方法再扩展一下。
先考虑朴素的情况。假如我们有两个 GF \(F(z),G(x)\),满足 \(F(z)\) 微分有限,它的 ODE 最短阶数为 \(n\)。我们现在就要求一个 \([x^{m-1}]F(G(x))\)。
首先,我们把复合构造成良复合形式。假设 \(c=[x^0]G(x)\),我们可以在 \(c\) 处展开 \(F(z)=\sum_{k\ge 0}f_k(z-c)^k\)。这样的话,在复合 \(G(x)\) 的时候我们就可以做一个截断 \(\mathscr F(z+c)=F(z+c)\bmod z^m\)。显然,该截断满足 \([x^{m-1}]\mathscr F(G(x))=[x^{m-1}]F(G(x))\)。展开之后,我们有 \([x^{m-1}]=[x^{m-1}]\sum_{k=0}^{m-1}([z^k] \mathscr F(z))G(x)^k\)。问题变成了如何快速计算 \(G(x)^k\) 的单项系数和 \(\mathscr F(z)\) 的所有系数。
接下来,我们可以用多种方式给出 \(\mathscr F\) 的系数。首先,如果 \(F\) 自己有好算的封闭形式,那么我们可以先使用多项式工业算出来 \(F(z+c)\bmod z^{m}\),然后做一个平移,这是朴素的用不着 ODE 的方法。
其次,我们用上之前对于截断的讨论。如果原先的 \(F\) 的 ODE 为:
我们可以通过某种方法找到一个修正因子 \(\mathscr D(z)\):
再平移回来得到:
现在就可以递推 \(\mathscr F\) 的系数了。
快速计算 \(G(x)^k\) 的单项系数的话......如果 \(G(x)\) 本身性质比较好,比如 \(G(x)=e^x\) 或者 \(\frac{1}{1-x}\) 之类的,我们可以很快地算出来单项系数。退一步,如果存在复合逆的话,我们也可以使用拉格朗日反演算 \([x^{m-1}]\frac{1}{1-yG(x)}\bmod y^m\)。其它情况的方法暂时不明。
实际上,我们可以将它推广到 \(G(x)^k\) 的线性组合,这也是 EI 在原博客中给出的形式:
如果对于任意的 \(0\le k<m\) 的整数 \(k\),下面的:
\[\sum_{j=0}^ma_j[x^j]G(x)^k \]都已知,那么我们可以通过上述方法快速地算出:
\[\sum_{j=0}^ma_j[x^j]F(G(x)) \]
从我们的结果到这一步是容易的,只需要交换求和号即可快速计算。
所以这个做法有一定的局限性:
-
由于我还不会什么 \(\mathscr D\) 的机械化算法,所以一旦 \(n\) 很大或者 \(p_k(z)\) 很长,这个方法就会变得尤其繁琐。
-
这个方法只能快速算一项结果。当需要一次求多项结果的时候,它会显得有点鸡肋,甚至可能直接失效......