Berlekamp–Massey 算法，Euclid 算法和有理函数重建

时隔多年把这个坑补一下。。如果有错误请告诉我！感谢。

设 $\deg 0=-\infty$。

有理函数重建

一说到这个我们就会想到寻找线性递推的 Berlekamp–Massey 算法，也有很多文章介绍了 Berlekamp–Massey 算法与 Euclid 算法是等价的，所以借助 Half-GCD 算法可以更快的完成这个求算。

$\mathbb{C}\left(\left( x^{-1}\right)\right)$ 上的乘法逆元

反形式 Laurent 级数环 $\mathbb{C}\left(\left( x^{-1}\right)\right)$ 可以被视为在多项式环 $\mathbb{C}\left\lbrack x\right\rbrack$ 上允许无限多项 $x$ 的指数为负数的项。类比于我们常见的 $\mathbb{Z}$ 和 $\mathbb{Q}$ 的关系一样。

这样一来，任一非零多项式在 $\mathbb{C}\left(\left( x^{-1}\right)\right)$ 我们都可以计算它的乘法逆元。因为 $\mathbb{C}$ 是域，所以 $\mathbb{C}\left(\left( x^{-1}\right)\right)$ 是域。

Euclid 除法是对 $A(x)$ 和非零多项式 $B(x)$ 计算

$$ A(x)=Q(x)B(x)+R(x) $$

其中 $\deg Q=\deg A-\deg B$ 且 $\deg R\lt \deg B$，那么 $Q(x)$ 就是 $A(x)$ 除以 $B(x)$ 的商数而 $R(x)$ 就是余数。

在环 $\mathbb{C}\left(\left( x^{-1}\right)\right)$ 上我们就可以写成

$$ \frac{A(x)}{B(x)}=Q(x)+\frac{R(x)}{B(x)} $$

而 $R(x)/B(x)$ 没有 $x^{\geq 0}$ 的项，且 $Q(x)$ 没有 $x^{\lt 0}$ 的项，这就是 Euclid 除法。

我们对于度数的定义有一个自然的扩展，即 $\deg\left(B^{-1}\right)=-\deg B$，同样的我们也可以定义向下取整函数 $\left\lfloor \cdot \right\rfloor :\mathbb{C}\left(\left( x^{-1}\right)\right)\to \mathbb{C}\left\lbrack x\right\rbrack$。

从反形式 Laurent 级数重建有理函数

考虑数列 $\left(a_j\right)_{j\geq 0}$ 是由有理函数 $P/Q$ 的系数生成的：$P/Q=\sum_{j\geq 1}a_{j-1}x^{-j}$，给出 $\left(a_j\right)_{j=0}^{k-1}$，求出 $P,Q$ 且最小化 $\deg Q$。我们发现

$$ \frac{\sum_{j=0}^{k-1}a_{k-1-j}x^j}{x^k}=\sum_{j=1}^{k}a_{j-1}x^{-j} $$

然后直接给出算法为（也许算法名不太准确，不过不要在意，后面懒得改了，似乎叫近似更好）

$$ \begin{array}{ll} &\textbf{Algorithm }\operatorname{RFR}(A,B,k)\text{:} \\ &\textbf{Input}\text{: }A,B\in\mathbb{C}\left\lbrack x\right\rbrack, \deg A\lt \deg B, k\in\mathbb{Z}_{\gt 0}\text{.} \\ &\textbf{Output}\text{: }C,D\in\mathbb{C}\left\lbrack x\right\rbrack\text{ such that } \\ &\qquad \qquad \left\lbrack x^{\left\lbrack -k,-1\right\rbrack}\right\rbrack A/B=\left\lbrack x^{\left\lbrack -k,-1\right\rbrack}\right\rbrack C/D, \\ &\qquad \qquad\deg C\lt \deg D,\deg D\text{ is minimized.} \\ 1&\textbf{if }\deg A-\deg B\lt -k\textbf{ then return }(0,1) \\ 2&Q\gets \left\lfloor B/A\right\rfloor,R\gets B\bmod{A} \\ &\text{Now we have }A/B=\cdots +c_{-k}x^{-k}+\cdots +c_{-\deg Q}x^{-\deg Q}\text{,} \\ &\text{If we set }A/B=1/(B/A)=1/(C+D)\text{ where} \\ &\qquad C,D\in\mathbb{C}\left(\left( x^{-1}\right)\right)\text{ and }\deg C=\deg Q\gt \deg D, \\ &\text{then we have }(A/B)\cdot C+(A/B)\cdot D=1\text{.} \\ &\text{If we have }\left\lbrack x^{\left\lbrack -k,-\deg Q\right\rbrack}\right\rbrack 1/C=\left\lbrack x^{\left\lbrack -k,-\deg Q\right\rbrack}\right\rbrack A/B\text{,} \\ &\text{then we can drop }D\text{. Which is saying that} \\ &\qquad (A/B)+(A/B)\cdot D/C=1/C \\ &\implies \deg A-\deg B+\deg D-\deg C\lt -k \\ &\implies \deg D\lt -k+2\deg Q \\ 3&\textbf{if }\deg R-\deg A\lt -k+2\deg Q \textbf{ then return }(1,Q) \\ &\text{If we set }C\gets Q\text{ and now }D=R/A\text{,} \\ &\deg D=\deg R-\deg A\lt -k+2\deg Q\text{,} \\ &\text{then we can just drop }D\text{.} \\ 4&(E,F)\gets \operatorname{RFR}(R,A,k-2\deg Q) \\ &\text{Otherwise we set }C\gets Q+E/F\text{ with }\deg F\text{ minimized.} \\ &\text{Now we have }\frac{1}{Q+E/F}=\frac{F}{QF+E}\text{.} \\ 5&\textbf{return }(F,QF+E) \end{array} $$

那么我们只需调用一次 $(A,B)\gets\operatorname{RFR}\left(\sum_{j=0}^{k-1}a_{k-1-j}x^j,x^k,k\right)$，那么 $A/B$ 就是我们想要的解。

接下来我们考虑算法的正确性，首先第一行和第三行返回时分母的度数一定是最小化的，只需考虑第五行返回时分母的度数是否最小化，我们意识到只要找到一个分母度数最小的有理函数生成了 $R/A$ 的 $x^{-k+2\deg Q},\dots,x^{-1}$ 的系数，那么简化之后的分母度数一定还是最小的，并且显然所有的返回值都是最简分数。

看完这个算法之后我们意识到，参数可以改为只有一个『项数有限的反形式 Laurent 级数』，也就是『有理函数截取一部分』也是可以的，但是这是否会使得计算量变大？

注意到第一步只有可能在第一次调用时返回而不可能在递归调用中返回，因为后面第三步没有返回时已经保证了 $R\neq 0$；在第三步我们真的需要返回 $\left(1,Q\right)$ 吗？如果不需要那么多项的话，$Q$ 的低位并不会对结果产生影响。但是使用这样的写法可以保证和 Euclid 算法几乎是一样的。

这个算法的时间和空间复杂度都是 $O\left(k^2\right)$，因为我们将 $R$ 放在栈上了，接下来考虑如何将算法改造为迭代。

如果我们将这个算法写下来，发现其是在计算连分数：

$$ \left\lbrack a_0,a_1,a_2,\dots\right\rbrack :=a_0+\cfrac{1}{a_1+\cfrac{1}{a_2+\cdots}} $$

设

$$ \frac{p_m}{q_m}=\left\lbrack a_0,\dots,a_m\right\rbrack $$

我们知道

$$ \frac{p_0}{q_0}=a_0,\quad \frac{p_1}{q_1}=\frac{a_0a_1+1}{a_1} $$

假设对于 $m\geq 2$ 有

$$ \frac{p_m}{q_m}=\frac{a_mp_{m-1}+p_{m-2}}{a_mq_{m-1}+q_{m-2}} $$

并且因为分子和分母互素，那么

$$ \begin{aligned} \left\lbrack a_0,\dots,a_m,a_{m+1}\right\rbrack&=\left\lbrack a_0,\dots,a_m+\frac{1}{a_{m+1}}\right\rbrack \\ &=\cfrac{\left(a_m+\cfrac{1}{a_{m+1}}\right)p_{m-1}+p_{m-2}}{\left(a_m+\cfrac{1}{a_{m+1}}\right)q_{m-1}+q_{m-2}} \\ &=\cfrac{a_mp_{m-1}+p_{m-2}+\cfrac{p_{m-1}}{a_{m+1}}}{a_mq_{m-1}+q_{m-2}+\cfrac{q_{m-1}}{a_{m+1}}} \\ &=\frac{a_{m+1}\left(a_mp_{m-1}+p_{m-2}\right)+p_{m-1}}{a_{m+1}\left(a_mq_{m-1}+q_{m-2}\right)+q_{m-1}} \\ &=\frac{a_{m+1}p_m+p_{m-1}}{a_{m+1}q_m+q_{m-1}} \\ &=\frac{p_{m+1}}{q_{m+1}} \end{aligned} $$

如果我们设 $p_{-1}=1,p_{-2}=0$ 和 $q_{-1}=0,q_{-2}=1$ 那么 $\frac{p_m}{q_m}=\frac{a_mp_{m-1}+p_{m-2}}{a_mq_{m-1}+q_{m-2}}$ 对于 $m\geq 0$ 都成立。

下面我们给出迭代的算法：

代码，有错误可以评论区指出吗！感谢！

和形式幂级数的关系

~~说到这里，相信聪明的你一定已经明白了！~~

对应的形式幂级数的有理函数就是 $\dfrac{x^{\deg Q-1}P\left(x^{-1}\right)}{x^{\deg Q}Q\left(x^{-1}\right)}$，注意此时分子的度数没有再保证小于分母的度数了。

利用 Half-GCD 算法求解

在 LOJ 写过一份 Yap 的 HGCD 算法的伪代码，根据定义，我们只需要调用一次 HGCD 算法，然后检查是否满足条件（观察我们上面的朴素算法每次 $k$ 会减少两倍的 $\deg Q$，HGCD 算法也正是执行“一半”的 GCD 算法，所以至多只需要一次调整就能够符合我们的要求），最多再进行一次 Euclid 除法就可以得到需要的有理函数了。

hly1204的博客

博客

有理函数重建

$\mathbb{C}\left(\left( x^{-1}\right)\right)$ 上的乘法逆元

从反形式 Laurent 级数重建有理函数

和形式幂级数的关系

利用 Half-GCD 算法求解

评论

发表评论