算法简介

Gosper 算法是一个用来求一些超几何项的离散不定积分的封闭形式的算法。特别地，所求的是一些形如 $$t_n=P(n)\cdot {\displaystyle \frac{\prod _{i=1}^s(a_{i}n+b_i)!}{\prod _{i=1}^t(c_{i}n+d_i)!}}$$ 的函数的离散不定积分。

形式化地，我们想要求出以另一个超几何项 $z_n$ 满足 $$\begin{array}{}\text{(1)}& z_{n+1}-z_n=t_n.\end{array}$$

本文大部分翻译自 Marko Petkovšek, Herbert S. Wilf 和 Doron Zeilberger 合著的 A=B 一书的 Chapter 5: Gosper's Algorithm。

记号

本文中 $\prod$ 和 $\sum$ 的优先级均大于 $\cdot$、$+$ 和 $\times$，小于隐式乘法。

在本文中，$\mathrm{deg}f$ 被定义为多项式 $f$ 的最高项次数，$\mathsf{lc}f$ 被定义为多项式 $f$ 的最高项系数，$\mathsf{sc}f$ 被定义为非常多项式 $f$ 的次高项系数。

本文用 $\perp$ 来描述“互素”，譬如 $f\perp g$ 等价于 $\mathrm{deg}gcd(f,g)=0$。

基本思路

首先我们不难发现 $${\displaystyle \frac{z_n}{t_n}}={\displaystyle \frac{z_n}{z_{n+1}-z_n}}={\displaystyle \frac{1}{\frac{z_{n+1}}{z_n}-1}}$$ 是一个关于 $n$ 的有理函数。于是我们设其为有理函数 $y(n)$，将 $z_n=y(n)t_n$ 代入 $(1)$ 式得到 $$y(n+1)t_{n+1}-y(n)t_n=t_n,$$ 设 ${\displaystyle \frac{t_{n+1}}{t_n}}=r(n)$ 则有 $$\begin{array}{}\text{(2)}& y(n+1)r(n)-y(n)=1.\end{array}$$ 假设我们能够把 $r(n)$ 写成 $$r(n)={\displaystyle \frac{a(n)c(n+1)}{b(n)c(n)}}$$ 的形式（在后文我们会给出一个算法来将任意有理函数 $r(n)$ 分解），满足 $$\begin{array}{}\text{(3)}& h\in \mathbb{N}⟹a(n)\perp b(n+h),\end{array}$$ 在后续的证明中，我们会发现我们将解有理函数的问题转为了解多项式的问题。

考虑将 $y(n)$ 写成 $$y(n)={\displaystyle \frac{b(n-1)x(n)}{c(n)}}$$ 的形式，代入 $(2)$ 式得 $$\begin{array}{}\text{(4)}& a(n)x(n+1)-b(n-1)x(n)=c(n).\end{array}$$ 显然 $x(n)$ 是有理函数。事实上我们能够证明，$x(n)$ 实际上甚至是多项式：

假设结论是错的，那么我们就能将 $x(n)$ 写成两个互素的多项式 $f(n)$ 和 $g(n)$ 之比 ${\displaystyle \frac{f(n)}{g(n)}}$，其中 $\mathrm{deg}g\ne 0$。代入 $(3)$： $$\begin{array}{}\text{(5)}& a(n)f(n+1)g(n)-b(n-1)f(n)g(n+1)=c(n)g(n)g(n+1).\end{array}$$ 设 $N$ 为使 $gcd(g(n),g(n+N))$ 不为常多项式的最大整数 $N$，容易知道 $N\ge 0$。设 $u(n)$ 为 $g(n)$ 和 $g(n+N)$ 的一个非常多项式的不可约公因式。将 $u(n-N)\mid g(n)$ 代入 $(4)$ 知： $$u(n-N)\mid b(n-1)f(n)g(n+1).$$ 因为 $u(n-N)\mid g(n)$，故 $u(n-N)\nmid f(n)$。若 $u(n-N)\mid g(n+1)$，则 $u(n)\mid g(n+N+1)$，即 $gcd(g(n+N+1),g(n))$ 不为常多项式，这与我们对 $N$ 的选取不符。因此 $u(n-N)\mid b(n-1)$，进而 $u(n+1)\mid b(n+N)$。

同样地，将 $u(n+1)\mid g(n+1)$ 代入 $(4)$ 式得： $$u(n+1)\mid a(n)f(n+1)g(n).$$ 由于 $u(n+1)\nmid f(n+1)$，以及 $u(n+1)\nmid g(n)$（因为 $u(n)$ 不能同时整除 $g(n-1)$ 和 $g(n+N)$），我们有 $u(n+1)\mid a(n)$。于是 $u(n+1)$ 是 $a(n)$ 和 $b(n+N)$ 的一个非常多项式公因子，这与假设 $(3)$ 矛盾，因而 $x(n)$ 是一个多项式。$◼$

现在我们确定了 Gosper 算法的基本流程：

1. 将 $r(n)={\displaystyle \frac{t_{n+1}}{t_n}}$ 分解成 ${\displaystyle \frac{a(n)c(n+1)}{b(n)c(n)}}$ 的形式，同时满足 $(3)$ 式的限制；
2. 找到 $(5)$ 式的解 $x(n)$，或发现无解；
3. 若有解，返回 $z_n={\displaystyle \frac{b(n-1)x(n)}{c(n)}}{t}_n$。

分解 $r(n)$

首先设 $r(n)={\displaystyle \frac{f(n)}{g(n)}}$（为了方便，令 $f$ 和 $g$ 均为首一多项式，这样可能会在最后需要乘一个系数），其中 $f\perp g$。如果 $a=f,b=g$ 满足 $(3)$ 的限制，那我们可以直接将 $r(n)$ 分解成 $a(n)=f(n)$，$b(n)=g(n)$，$c(n)=1$ 的形式。

否则，如果存在 $h\in \mathbb{N}$ 使得 $f(n)$ 和 $g(n+h)$ 有公因式 $u(n)$，则我们可以从 $f(n)$ 中提取一个因式 $u(n)$，从 $g(n)$ 中提取一个因式 $u(n-h)$，从而递归下去。分解出来的因式直接让 $c(n)$ 乘上 $\prod _{j=1}^{h}u(n-j)$ 即可。

唯一问题是，如何判断 $f$ 和 $g$​ 是否满足 $(3)$ 式？

一个做法是使用多项式结式（resultant）。

域 $F$ 上两个关于 $n$ 的多项式 $A(n)$ 和 $B(n)$ 的结式 $\mathrm{res}(A,B)$ 由下式定义： $$\mathrm{res}(A,B):=(\mathsf{lc}B{)}^{\mathrm{deg}A}(\mathsf{lc}A{)}^{\mathrm{deg}B}\prod _{(\mu ,\lambda )\in {\bar{F}}^2\wedge A(\lambda )=0\wedge B(\mu )=0}(\mu -\lambda ).$$ 其中域 $\bar{F}$ 是域 $F$ 的代数闭包，即，所有以 $F$ 中元素为系数的多项式的解都在 $\overline{F}$ 内。

事实上，$\mathrm{res}(A,B)$ 的一个等价定义是： $$\mathrm{res}(A,B)=(\mathsf{lc}B{)}^{\mathrm{deg}A}\prod _{\mu \in \bar{F}\wedge B(\mu )=0}A(\mu ).$$ 根据以上两种定义不难得到结式的以下性质：

• $\mathrm{res}(A,B)=(-1{)}^{\mathrm{deg}A\cdot \mathrm{deg}B}\prod \mathrm{res}(B,A)$；
• $\mathrm{res}(\mu A,\lambda B)={\mu }^{\mathrm{deg}A}{\lambda }^{\mathrm{deg}B}\mathrm{res}(A,B)$；
• $\mathrm{res}(A,B)=(\mathsf{lc}B{)}^{\mathrm{deg}A}(\mathsf{lc}A{)}^{\mathrm{deg}B}$，若 $\mathrm{deg}A\cdot \mathrm{deg}B=0$；
• $\mathrm{res}(A-CB,B)=\mathrm{res}(A,B)$，若 $\mathsf{lc}B=1$。

于是我们可以得到一个 $O(\mathrm{deg}A\cdot \mathrm{deg}B)$ 的求解结式的 Euclid 算法：

• 若 $\mathrm{deg}B=0$，则 $\mathrm{res}(A,B)=(\mathsf{lc}B{)}^{\mathrm{deg}A}(\mathsf{lc}A{)}^{\mathrm{deg}B}$；
• 否则，$\mathrm{res}(A,B)=(\mathsf{lc}B{)}^{\mathrm{deg}A-\mathrm{deg}(AmodB)}(-1{)}^{\mathrm{deg}B\cdot \mathrm{deg}(AmodB)}\mathrm{res}(B,AmodB)$。

事实上，还有一种被称为 Half-GCD 的算法可以在 $O(n^{1+\epsilon })$ 的时间复杂度内求出 $A$ 和 $B$ 的结式，其中 $n=max\{A,B\}$。${}^{[1]}$

不难发现，两个非常多项式 $f$ 和 $g$ 有公共根当且仅当 $\mathrm{res}(f,g)=0$。我们设 $R(h)=\mathrm{res}(f(n),g(n+h))$，则所有违反 $(3)$ 式的 $h$ 即为 $R(h)$ 的非负整数根。

如何在 $\mathbb{N}$ 上解 $R(h)=0$？首先我们仅保留和分母互素的分子 $S(h)$，于是只需要解 $S(h)=0$，其中 $S(h)$ 是多项式。如果 $S$ 的系数在 $\mathbb{Q}$ 内，则我们可以通过先将其化为首一多项式，并枚举其常数项的因数来得到其所有的有理根。事实上，还有一种更加高效的、利用 $p$ 进数的算法可以解决这个问题。${}^{[2]}$

更加通用地，设 $K$ 为以 $\bar{K}$ 中元素为系数的关于 $\alpha$ 的有理函数所组成的元素。假设我们已经拥有一种处理 $S$ 的系数在 $A_{\bar{K}}$ 内的算法，那么我们可以以如下方式构建出 $A_K$ 所对应的算法：

首先做转换 $$\begin{array}{rl}R(h)& =\sum _{i=0}^s{a}_i{h}^i\\ & ={\displaystyle \frac{1}{r(\alpha )}}\sum _{i=0}^s{p}_i(\alpha )h^i\\ & ={\displaystyle \frac{1}{r(\alpha )}}\sum _{i=0}^s{h}^i\sum _{j=0}^{t_s}{\mu }_{i,j}{\alpha }^j\\ & ={\displaystyle \frac{1}{r(\alpha )}}\sum _{j=0}^t{\alpha }^j\sum _{i=0}^s{\mu }_{i,j}{h}^i\\ & ={\displaystyle \frac{1}{r(\alpha )}}\sum _{j=0}^t{q}_j(h){\alpha }^j,\end{array}$$ 其中 $r,p_i,q_j\in \bar{K}[x]$。如果有 $u\in \bar{K}$ 满足 $R(u)=0$，那么就有 $$\begin{array}{}\text{(6)}& \sum _{j=0}^t{q}_j(u){\alpha }^j=0.\end{array}$$ 分类讨论：

• 如果 $\alpha$ 在 $\bar{K}$ 中是超越的，那么显然 $q_j(u)=0$ 对全部 $[0,t]$ 内的整数 $j$ 成立；
• 如果存在一个系数在 $\bar{K}$ 内的多项式的根为 $\alpha$，且 $d$ 为多项式的最小次数。则 $\mathrm{deg}{p}_i<d$。于是，$t<d$。于是使 $(6)$ 式成立的唯一可能是 $q_j(u)=0$ 对全体 $[0,t]$ 内整数 $j$ 成立。

无论如何都有 $q_j(u)=0$。于是我们可以对 $q_j$ 分别应用 $A_{\bar{K}}$。

另一种做法是将 $f$ 和 $g$ 化成若干个不可约多项式乘积（注意不必在整个代数闭包上完全分解），那么我们需要找到一个 $f(n)$ 的不可约因式 $u(n)$ 和 $g(n)$ 的一个不可约因式 $v(n)$ 满足存在 $h\in \mathbb{N}$ 使得 $u(n)=v(n+h)$，这样 $f(n)$ 和 $g(n+h)$ 就有一个非常多项式公因式，即 $h$ 违反 $(3)$ 式。注意到一个因式对 $(u,v)$ 满足条件的一个必要条件是 $\mathrm{deg}u=\mathrm{deg}v=d$。设 $u(n)=n^d+A{n}^{d-1}+O(n^{d-2})$，$v(n)=n^d+B{n}^{d-1}+O(n^{d-2})$，那么通过比较次大项系数我们发现唯一可能满足条件的 $h={\displaystyle \frac{A-B}{d}}$。

在实践中，通常 $f$ 和 $g$ 已经被分解成一次因式了，这可能是因为初始的超几何项即为二项式系数和阶乘的乘积。此时，两种做法是完全相同的。如果 $f$ 和 $g$ 未被分解，则通常而言就算要使用结式，事先分解 $f$ 和 $g$ 也通常能加速结式的计算。那为什么还需要使用结式做法呢？这是因为在任何域上我们都能计算多项式结式，但目前人们还尚未知道通用的因式分解方法。

不难验证，通过“基本思路”中过程得到的 $r(n)$ 的分解 ${\displaystyle \frac{a(n)c(n+1)}{b(n)c(n)}}$ 是一个互素、唯一、次数最小的形式。

解 $x(n)$

在这节我们希望能解出满足 $(4)$ 式的多项式 $x(n)$。假定我们能知道 $\mathrm{deg}x(n)$ 的一个范围，我们就可以使用待定系数法列出一个线性方程组以解出 $x(n)$。

设 $\mathrm{deg}x(n)=d$，我们希望知道 $d$ 的可能值。可以分类讨论一下：

• 如果 $\mathrm{deg}a\ne \mathrm{deg}b$ 或 $\mathsf{lc}a\ne \mathsf{lc}b$，即 $a$ 和 $b$ 的最高项不同，那么 $(4)$ 式左手边的系数为 $d+max\{\mathrm{deg}a,\mathrm{deg}b\}$，右手边的系数为 $\mathrm{deg}c$。$d$ 的唯一可能值为 $\mathrm{deg}c-max\{\mathrm{deg}a,\mathrm{deg}b\}$；

• 否则，即 $a$ 和 $b$ 的最高项相同。即 $(4)$ 式左手边的最高项被消去了。继续分类：

  • 如果 $(4)$ 式左手边的次高项未被消去，则 $(4)$ 式左手边的系数为 $d+max\{\mathrm{deg}a,\mathrm{deg}b\}+1$，右手边的系数为 $\mathrm{deg}c$。$d$ 的唯一可能值为 $\mathrm{deg}c-max\{\mathrm{deg}a,\mathrm{deg}b\}-1$；

  • 否则，设 $$\begin{array}{rl}a(n)& =\lambda n^k+A{n}^{k-1}+O(n^{k-2}),\\ b(n-1)& =\lambda n^k+B{n}^{k-1}+O(n^{k-2}),\\ x(n)& =C_0{n}^d+C_1{n}^{d-1}+O(n^{d-2}).\end{array}$$ 代入 $(4)$ 式得到： $$a(n)x(n+1)-b(n-1)x(n)=C_0(\lambda d+A-B)n^{k+d-1}+O(n^{k+d-2}).$$ 于是 $d$ 的唯一可能值为 ${\displaystyle \frac{B-A}{\lambda }}$。

  于是，此情况中 $d$ 仅可能在 $\{\mathrm{deg}c-max\{\mathrm{deg}a,\mathrm{deg}b\}-1,\frac{\mathsf{sc}(b(n-1))-\mathsf{sc}a}{\mathsf{lc}a}\}$ 中取值。

先得到 $d$ 的可能取值范围，找到其中最大的自然数 $d_m$。如果 $d_m$ 不存在则无法找到封闭形式，否则直接待定系数法。如果待定系数法无解则亦无法找到封闭形式。

小练习

通过 Gosper 算法求下述和式的封闭形式。 $$\begin{array}{ll}\text{(a)}& {\displaystyle \sum _{n=0}^m{\displaystyle \frac{{(\genfrac{}{}{0ex}{}{2n}{n})}^2}{4^{2n}(n+1)}}}\\ \text{(b)}& {\displaystyle \sum _{n=0}^m{\displaystyle \frac{(4n-1){(\genfrac{}{}{0ex}{}{2n}{n})}^2}{4^{2n}(2n-1{)}^2}}}\\ \text{(c)}& {\displaystyle \sum _{n=0}^m{\displaystyle \frac{n(n-\frac{1}{2})!}{(n+1)!}}}\\ \text{(d)}& {\displaystyle \sum _{n=0}^m{\displaystyle \frac{n(n+a+b)a^n{b}^n}{(n+a)!(n+b)!}}}\\ \text{(e)}& {\displaystyle \sum _{n=1}^m{\displaystyle \frac{6n+3}{4n^4+8n^3+8n^2+4n+3}}}\\ \text{(f)}& {\displaystyle \sum _{n=4}^m{\displaystyle \frac{4(1-n)(n^2-2n-1)}{n^2(n+1{)}^2(n-2{)}^2(n-3{)}^2}}}\\ \text{(g)}& {\displaystyle \sum _{n=1}^m{\displaystyle \frac{2^n(n^4-14n^2-24n-9)}{n^2(n+1{)}^2(n+2{)}^2(n+3{)}^2}}}\\ \text{(h)}& {\displaystyle \sum _{n=0}^m{\displaystyle \frac{n^4{4}^n}{(\genfrac{}{}{0ex}{}{2n}{n})}}}\end{array}$$

其中 $(\mathrm{a}),(\mathrm{b}),(\mathrm{c}),(\mathrm{d}),(\mathrm{h})$ 摘自 ${}^{[3]}$，$(\mathrm{e})$ 摘自 ${}^{[4]}$，$(\mathrm{f}),(\mathrm{g})$ 摘自 ${}^{[5]}$。

参考资料

• Marko Petkovšek, Herbert S. Wilf and Doron Zeilberger, A=B, 83-100.

致谢

感谢 @Y_B_X 同学给予的许多指导和帮助。

Reference

${}^{[1]}$: E.I., HALF-GCD 算法的阐述.

${}^{[2]}$: Loos, R., Computing rational zeros of integral polynomials by $p$-adic expansion, SIAM J. Comp. 12 (1983), 286-293.

${}^{[3]}$: Gosper, R. W., Jr., Indefinite hypergeometric sums in MACSYMA, in: Proc. 1977 MACSYMA Users' Conference, Berkeley, 1977, 237-251.

${}^{[4]}$: Abramov, S. A., On the summation of rational functions, USSR Comp. Maths. Math. Phys. 11 (1971), 324-330.

${}^{[5]}$: Man, Y. K., On computing closed forms for indefinite summations, J. Symbolic Computation 16 (1993), 355-376.