Berlekamp–Masseyのアルゴリズム

問題：形式的べき級数 $F(x) = \sum_{i=0}^\infty a_i x^i$ が分数 $A(x) / B(x)$ の形で書けるとする。$\deg A < n, \deg B \leq n$ が分かっている。$F$ の先頭 $2n$ 項を切り取った $f(x) := F(x) \bmod x^{2n}$ だけが与えられる。$A(x), B(x)$ を求めよ。

$A, B$ を定数倍して $B$ を monic にできます。従って、未知数は $2n$ 個あります。分母を払って $f(x)B(x) = A(x) \pmod {x^{2n}}$ の両辺の係数を比較すれば $2n$ 個の線形方程式が得られ、ガウスの消去法 $O(n^3)$ で $A, B$ が得られます。

Berlekamp–Masseyのアルゴリズムでは $O(n^2)$ 時間にできます。Half-GCD を使うと更に $O(n(\log n)^2)$ にできます。

このような有理多項式は、存在するなら、一意に定まります。$f_1 / g_1 = f_2 / g_2$ より $f_1 g_2 = f_2 g_1$ で、次数は $2n$ 未満です。$f_1$ と $g_1$、$f_2$ と $g_2$ が互いに素、$g_1,g_2$ が monic であることより、$f_1 = f_2, g_1 = g_2$ です。

※ 解が存在しない場合があります。例えば、$\frac{a}{1+bx} = x \pmod {x^2}$ の分母を払って各係数を比較すると $a = x$ ですが、$a$ は定数なので矛盾。

1 参考文献
2 連分数展開
- 2.1 整理できていないこと
3 ユークリッドの互除法
4 $\mathbb{F}_p$ の元の分数表示

連分数を使ったアルゴリズムと、ユークリッドの互除法を使ったアルゴリズムの二つを紹介します。

参考文献

下記文献を参考にしました。[1] は私が最初に Berlekamp-Massey のアルゴリズムを勉強した文献です。オリジナルの方法とは少々違っていて、ユークリッドの互除法を全面に使っています。[3] はBerlekamp-Massey のアルゴリズムを連分数から見た解説です。[3] の細かい理論を確認するために [2] を読みました。

[1] Sugiyama, Yasuo, et al. A method for solving key equation for decoding Goppa codes. Information and Control 27.1 (1975): 87-99.
[2] Mills, William H. “Continued fractions and linear recurrences.” Mathematics of Computation 29.129 (1975): 173-180.
[3] かならいさんの連分数の観点からの解説

連分数展開

$F(x)$ を連分数展開することにより、$A, B$ を $O(n^2)$ で求めることができます。

かならいさんの解説が素晴らしいので、そちらも見てください。

整数部分を引きはがし、小数部分の逆数を取ることを繰り返した形を連分数展開と言います。例えば、円周率 $\pi$ のとき、次のようになります。
$$\begin{align}
\pi & = 3.14159265358979\ldots \\
& = 3 + 0.14159265358979\ldots\\
& = 3 + \frac{1}{\frac{1}{0.14159265358979\ldots}} \\
& = 3 + \frac{1}{7.06251330593104\ldots} \\
& = 3 + \frac{1}{7+0.06251330593104\ldots} \\
& = 3 + \frac{1}{7+\frac{1}{\frac{1}{0.06251330593104\ldots}}} \\
& = 3 + \frac{1}{7+\frac{1}{15.9965944066857\ldots}} \\
& = 3 + \frac{1}{7+\frac{1}{15+0.9965944066857\ldots}}
\end{align}$$
$0.9965944066857\cdots\approx 1$ として連分数展開を途中で打ち切った $3 + \frac{1}{7+\frac{1}{16}} = \frac{355}{113} = 3.1415929203\ldots$ は $\pi$ の非常に良い近似になっています。実は、連分数展開を途中で打ち切ってごにょごにょした値は最良近似であることが知られています。

有理数 $\alpha$ と自然数 $n$ が与えられる。分母が $n$ 以下の下で $|\alpha-\beta|$ を最小化する $\beta$ が連分数展開により $O(\log n)$ で求められる (Khintchine [1956])。

証明は組合せ最適化第2版 (理論とアルゴリズム)の定理4.8をご参照ください。

$x$ の冪が大きいほどノルムが小さいと考えて、有理数の場合と同様に $xF$ を連分数展開します。
$$\begin{align}
xF(x) &= x+x^2+2x^3+3x^4+5x^5+\cdots\\
&= \frac{1}{\frac{1}{x+x^2+2x^3+3x^4+5x^5+\cdots}} \\
&= \frac{1}{x^{-1}-1+\frac{-x^2-x^3-2x^4-3x^5-5x^6-\cdots}{x+x^2+2x^3+3x^4+5x^5+\cdots}} \\
&= \frac{1}{x^{-1}-1+\frac{1}{\frac{x+x^2+2x^3+3x^4+5x^5+\cdots}{-x^2-x^3-2x^4-3x^5-5x^6-\cdots}}} \\
&= \frac{1}{x^{-1}-1+\frac{1}{-x^{-1}}} \\
\end{align}$$

$xF$ の近似値である $xf$ を連分数展開します。
$$\begin{align}
xf(x) &= x+x^2+2x^3+3x^4 \\
&= \frac{1}{\frac{1}{x+x^2+2x^3+3x^4}} \\
&= \frac{1}{x^{-1}-1+\frac{-x^2-x^3+3x^4}{x+x^2+2x^3+3x^4}} \\
&= \frac{1}{x^{-1}-1+\frac{1}{\frac{x+x^2+2x^3+3x^4}{-x^2-x^3+3x^4}}} \\
&= \frac{1}{x^{-1}-1+\frac{1}{-x^{-1}+\frac{5x^3+3x^4}{-x^2-x^3+3x^4}}} \\
\end{align}$$
$xf$ の連分数展開は $xF$ の連分数展開と誤差項 $\frac{5x^3+3x^4}{-x^2-x^3+3x^4}$ で表せました。誤差項を消せば、$xF$ の連分数展開になります。$\bmod {x^{2n}}$ で打ち切ったとき、常にこのようなことが成り立ちます。ちなみに、$F$ と $f$ の連分数展開は一致するとは限りません。

$$\begin{align} F &= 1 + \frac{1}{x^{-1}-2+\frac{1}{x^{-1}+1}} \\ f& = 1 + \frac{1}{x^{-1}-2+\frac{1}{x^{-1}-4+\frac{27x^3}{x^2+6x^3}}} \end{align}$$
連分数展開が得られれば、分数を末尾から解消することで全体を一つの分数で表せます（実際のアルゴリズムでは $q_i$ が得られるたびに、漸化式に従って逐次更新します）。では、連分数展開が実際に一致することを示しましょう。

定理：$xf$ の連分数展開は、$xF$ の連分数展開を完全に含む。

$f$ を $x^{-1}$ の多項式としてみたときの次数を $\deg^{-1} f$ と書くことにします。

$xF$ が次のように連分数展開されたとします。
$$ xF = \frac{1}{q_0+\frac{1}{q_1+\frac{1}{\ddots+\frac{1}{q_m} }}}$$
ただし、$$\begin{align}
r_0 &= 1 \\
r_1 &= xF \\
r_{i+2} &= r_{i}-q_{i} r_{i+1} \ (\deg^{-1} r_{i+2} < \deg^{-1} r_i) \\
\end{align}$$
とします。$\pi \approx \frac{355}{133}$ を求めたときと同様に、連分数展開の末尾の項から順に分数を解消することで、分母と分子が互いに素な分数が得られます（行列で書くとユニモジュラになることから従います。証明略）。

$d_i := \deg^{-1} q_i =\deg^{-1} r_i-\deg^{-1} r_{i+1}$ です。$\Delta_i := \sum_{k=0}^i d_k$ と置きます。$x^{-1}$ の多項式としてみたときの分母分子の次数は $\Delta_m,\Delta_m-d_0$ です。解が存在するなら、分母には定数項が存在します（証明は[2]の定理２）。従って $\Delta_m \leq n$ です。$d_0=0,\Delta_m=n$ のとき分子の $x$ の次数が $n$ になってしまいますが、求めているのは $xF$ なので問題ありません。

$q_i$ を計算するのに必要なのは、$x^{-1}$ の多項式としてみたときの、$r_i, r_{i+1}$ の上位 $1+d_i$ 項だけです（これは Half-GCD の中核的アイディアでした）。$xf$ での連分数展開において、$r_i$ の上位何桁が完全に $xF$ のそれと一致する形で残っているでしょうか。最初に $x$ を $f$ に掛けたので、項は高々 $x^{0},x^1,\ldots,x^{2n}$ まであります。ちょうど上から $x^{0},x^1,\ldots,x^{\Delta_{i-1}-1}$ は $0$ になっています。また、下から高々 $x^{2n},x^{2n-1},\ldots,x^{2n+1-\Delta_{i-2}}$ については、各 $k$ について $q_{k} r_{k+1}$ を引く際に、$r_{k+1}$ の下位の桁が $\bmod {x^{2n+1}}$ で削られて足りないがために、正確な値になっていない可能性があります。従って、確実に $xf$ と $xF$ で一致するのは $r_i$ の上から
$$\begin{align} &\ 2n+1-\Delta_i-\Delta_{i-1} \\
= &\ 2n+1-(\Delta_{i-1}+d_i)-\Delta_{i-1} \\
= &\ 2n-2\Delta_{i-1} + d_i+1 \\
\geq &\ d_i+1
\end{align}$$項だけです。これは、$q_i$ を完全な形で計算するのに十分。従って、$xf$ は $xF$ の連分数展開を完全な形で含みます。

整理できていないこと

$xF$ の連分数展開を途中で打ち切って、$x^{-1}$ の有理式として整理した形を $\frac{a_i}{b_i} = \frac{1}{q_0+\frac{1}{q_1+\frac{1}{\ddots+\frac{1}{q_i} }}}$ と置きます。帰納法により
$\begin{pmatrix}
q_i & 1 \\
1 & 0
\end{pmatrix}\cdots\begin{pmatrix}
q_1 & 1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix}
q_0 & 1 \\
1 & 0
\end{pmatrix}\begin{pmatrix}
a_{-1} & b_{-1} \\
a_{-2} & b_{-2}
\end{pmatrix}=\begin{pmatrix}
a_{i+1} & b_{i+1} \\
a_i & b_i
\end{pmatrix}$ の関係が成り立ちます。
$a_{-2}/b_{-2}=1/0, a_{-1}/b_{-1}=0/1$ なので、右辺の最後の行列は単位行列になります。両辺の行列式を比べることで、$a_{i+1}b_i-b_{i+1}a_i=(-1)^i$ を得ます。よって$$\frac{a_{i+1}}{b_{i+1}}-\frac{a_{i}}{b_{i}}= \frac{(-1)^{i}}{b_ib_{i+1}} $$となります。$1/b_ib_{i+1}$ を $x$ の有理式として整理すると、分子が $x^{\deg_{-1}(b_i)+\deg_{-1}(b_{i+1})}b_ib_{i+1}$ となります。従って、$\deg_{-1}(b_i)+\deg_{-1}(b_{i+1}) \geq n+1$ となる最小の $i$ を用いると $$xF = \frac{a_i}{b_i} \pmod {x^{n+1}}$$ となります。$ a_i / b_i $ を $x$ の有理式として整理すると、$\deg a_i \geq \deg b_i$ かつ $b_i$ の定数項が非零である必要がありますが、常にそうなるようです（証明はよすぽジャッジ……）。

ユークリッドの互除法

互除法を全面に使った方法でも解くことができます。$f(x)B(x) = A(x) \bmod x^{2n},\deg A \leq n, \deg B \leq n-1$ を満たす $A, B$ が求まればいいわけですが、$$B(x)f(x) + C(x)x^{2n} = A(x)$$の形にすれば、$f(x)$ と $x^{2n}$ の互除法の形になっています。互除法の $i$ 回目の iteration では、$B_i f + C_i x^{2n} = A_i$ だとします。

互除法を途中で打ち切ることにより、$\deg A_i < n-1, \deg B_i \leq n$ となる瞬間が必ず存在します。これは $\deg A_{i-1} \geq n, \deg A_i < n$ のとき
$$\begin{align} \deg B_i &\leq \sum_{k=0}^i (\deg A_{k-1}-\deg A_k) \\ &=2n-\deg A_{i-1} \\ &\leq n\end{align}$$であることから分かります。よって互除法を適切なタイミングで打ち切ることで、$O(n^2)$ 時間で $A, B, C$ が求まります。さらに、この $A, B$ は、有理式の表示が存在するなら、互いに素になっています（証明）。$B$が非可逆の場合、$\frac{A}{B}$ の形で書くことはできません。

$\mathbb{F}_p$ の元の分数表示

noshi91さんが $\mathbb{F}_p$ の元に対しても同じようなことができると言っていたので、考えてみました。

素数 $p$ と $a \in \mathbb{F}_p$ が与えられる。$\gcd(x,y)=1,|x| \leq \sqrt p, y < \sqrt p$ の下、$$\frac{y}{x} = a \pmod p$$ を満たす $x,y$ を求めよ。

残念ながら、$-8 / 5 = 7 / 9 = 84 \pmod {107}$ のように一意な表現にはなりません（$a \neq 0$ に対して、表現は高々二通りしか存在しません）。

$(x_0,z_0)=(0,1),(x_1,z_1)=(1,0)$ として $ax_i + pz_i = y_i$ で互除法をします。$y_{n} \geq \sqrt p, y_{n+1} < \sqrt p$ となる $n$ が存在します。帰納法（by noshi91さん）より $x_n, x_{n-1}$ が異符号なので $y_nq_n \leq y_{n-1}$ を用いると、$$\begin{align}
&y_n(-q_nx_{n}+x_{n-1})=y_nx_{n+1}\\
\Rightarrow &|-y_{n-1}x_{n}+y_nx_{n-1}| \geq y_n|x_{n+1}|
\end{align}$$です。$\begin{pmatrix}
q_i & 1 \\
1 & 0
\end{pmatrix}\cdots\begin{pmatrix}
q_2 & 1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix}
q_1 & 1 \\
1 & 0
\end{pmatrix}\begin{pmatrix}
x_{1} & y_{1} \\
x_{0} & y_{0}
\end{pmatrix}=\begin{pmatrix}
x_{i+1} & y_{i+1} \\
x_i & y_i
\end{pmatrix}$
なので $\det$ を取ることで、$x_{i+1}y_i-y_{i+1}x_i = (-1)^i(x_1y_0-y_1x_0) = (-1)^ip$ です。従って
$$p \geq y_n|x_{n+1}|$$を得ます。$x_{n+1} \leq p / y_n \leq \sqrt{p}, y_{n+1} < \sqrt{p}$ なので条件を満たす $x, y$ が得られました。また、以上より任意の $a \in \mathbb{F}_p$ に対して、このような分数表示が可能であることが分かりました。

code

連分数展開の方でも同じようにできるのでしょうか？