数列の漸化式の特定と第N項の計算
以下の数列があるとする。
- 最初の何項か $a_1,a_2,...,a_M$ がわかっている
- あんまり多くない次数の線形漸化式で表現できることがわかっている（具体的な式は分からなくてもいい）
- ○線形 $a_n = 3a_{n-1} + 5a_{n-2} + 7a_{n-3}$（この場合は $3$ 次）
- ×線形でない $a_n = a_{n-1}^2$
- ×線形でない $a_n = n \times a_{n-1}$
この時、以下のことが出来る。
- Berlekamp-Massey（バーレカンプ-マッシー）のアルゴリズム
- 特性方程式を特定する
- 漸化式を $a_n = 3a_{n-1} + 5a_{n-2}$ とすると、$x^2 - 3x - 5 = 0$ みたいになるやつ
- ここから漸化式の形もすぐ求められる
- あり得るもののうち、次数が最小のものが求まる
- そのため、はじめに与える数列の長さが足りないと、本来の特性方程式が求まる保証は無い
- 正解が $d$ 次漸化式の場合、$M=2d$ 項あれば十分
- Bostan-Mori のアルゴリズム
- 母関数 $P(x)/Q(x)$ で表される数列の第 $N$ 項を高速に計算する
- 母関数は特性方程式から導出できる
- $P,Q$ を $d$ 次多項式として、$O(d^2 \log{N})$ または $O(d \log{d} \log{N})$
この2つを続けて使用することで、最初の何項かがわかっているだけの数列の、かなり大きな $N$ に対する第 $N$ 項を求められる。
Berlekamp-Massey algorithm
イメージとしては、ユークリッドの互除法を、数列（形式的べき級数）に拡張して適用したもの。
- 広島大学 先端数学講義資料 6p.～
入力
数列の最初の何項か $A=(a_1,a_2,...,a_M)$
出力
$A$ の特性方程式。
より正確には、長さ $L+1$ の数列 $C$ であって、全ての $L \le i \le M$ に対して以下を満たすもの。
- $C_0a_i+C_1a_{i-1}+C_2a_{i-2}+...+C_La_{i-L}=0$
ただし、$C_0=1$ とする。（全ての項に同じ数を掛けたものを考えると答えが一意に定まらないため）
擬似コード
英語版Wikiに載ってる。
# Python寄りの擬似コード
# 初期化
C = [1] # 求める数列
B = [1] # 1つ前のCの状態を保存
L = 0 # Cの長さ-1
m = 1 # ポインタ？っぽいもの
b = 1 # 前回のdの値
for n in range(len(A)):
d = C[0]*A[n] + C[1]*A[n-1] + ... + C[L]*A[n-L]
# ①d=0なら、現在のCで、An を求める漸化式は成り立っていることになる
if d == 0:
m += 1
# それ以外なら調整する
# ②Cを拡張する場合
elif 2 * L <= n:
T = C
C -= d/b * (Bをmだけ右シフトしたもの)
L = n + 1 - L
B = T
b = d
m = 1
# ③拡張しない場合
else:
C -= d/b * (Bをmだけ右シフトしたもの)
m += 1
# これで C が特性方程式の係数になっている
$C$ が拡張される（$L$ が増加する）と $A$ のうち所与として扱われる範囲が変わるので、 多分その辺の関係で②と③で $B,b,m$ の更新の有無が変わってるんだと思う（イマイチよくわかってない）。
Bostan-Mori algorithm
前準備
Berlekamp-Massey の出力は特性方程式だったが、Bostan-Mori への入力は母関数 $\dfrac{P(x)}{Q(x)}$ の $P,Q$（の係数配列）でなくてはならない。
- Berlekamp-Massey の出力 $C$ は、そのまま $Q$ となる。
- 「$Q$」と、Berlekamp-Massey の入力でもあった「$A$ の最初の $d$ 項」を畳み込み、最初の $d$ 項を取ると $P$ となる。
たとえば
A[1] = A[2] = A[3] = 1
A[n] = A[n-1] + 2*A[n-2] + 3*A[n-3] のとき
Berlekamp-Massey の結果
Q: [1, -1, -2, -3]
Qと最初の3項を（たたみ込みの容量で）かけ、最初の3項をPとする
Q×[1, 1, 1]: [1, 0, -2, -6, -5, -3]
↓
P: [1, 0, -2]
すると、$\dfrac{P}{Q}$ は $A$ の母関数となっている
__1___1___1___6__11___... ← ここにAが現れる
[1,-1,-2,-3] ) 1 0 -2
1 -1 -2 -3
1 0 3 0
1 -1 -2 -3
1 5 3 0
1 -1 -2 -3
6 5 3 0
6 -6 -12 -18
11 15 18 0
11 -11 -22 -33
27 40 33 0
...
入力
- $A$ の母関数を表す、高々 $d$ 次の多項式 $P,Q$
- 求めたい項 $N$（0-indexed）
出力
- 数列 $A$ の第 $N$ 項
擬似コード
説明中の数列の「○○項目」は0-indexとする。
P = [1, 0, -2]
Q = [1, -1, -2, -3]
while N:
Q1 = (Qの奇数項目の正負を逆転した数列)
# たとえばループ1回目なら Q1 = [1, 1, -2, 3]
PQ = P * Q1 # たたみ込み
QQ = Q * Q1 # たたみ込み
if N & 1 == 0:
P = (PQの偶数項目のみを取りだした数列)
else:
P = (PQの奇数項目のみを取りだした数列)
Q = (QQの偶数項目のみを取り出した数列)
N >>= 1
# この時点の P[0]/Q[0] が答え
$Q$ として Berlekamp-Massey の結果を使った場合、初期の $Q_0=1$ となっている。 計算の過程上、これは不変なので、最後の答えの出力では、単に $P_0$ を参照するだけでもよい。
たたみ込みを愚直にやると $O(d^2 \log{N})$、高速フーリエ変換などを使えるなら $O(d\log{d}\log{N})$ となる。
やっていることのイメージ
$\dfrac{P(x)}{Q(x)}$ の $x^N$ の係数を求めたい。
分子と分母に $Q(-x)$ をかける。（上記擬似コードの Q1 を畳み込む操作）
- $\dfrac{P(x)Q(-x)}{Q(x)Q(-x)}$
すると、$Q(x)Q(-x)$ は、$x^1,x^3,x^5,...$ の項は全て打ち消され、$x^0,x^2,x^4,...$ についての項のみ残る。
分母が $x^{偶数}$ の項しかないので、
- $\dfrac{P(x)}{Q(x)}$ の $x^{偶数}$ の項には、$P(x)Q(-x)$ の $x^{偶数}$ の項のみが影響し、奇数の項は影響しない。
- $\dfrac{P(x)}{Q(x)}$ の $x^{奇数}$ の項には、$P(x)Q(-x)$ の $x^{奇数}$ の項のみが影響し、偶数の項は影響しない。
よって、以下のようにすると
- $U(x) = Q(x)Q(-x)$ の $x$ の指数を1/2にした多項式
- $V(x) = P(x)Q(-x)$ の $x^{偶数}$ の項のみ抽出し、指数を1/2にした多項式
- $W(x) = P(x)Q(-x)$ の $x^{奇数}$ の項のみ抽出し、指数から1引いて1/2にした多項式
このように変形できる。
- $\dfrac{P(x)}{Q(x)}=\dfrac{V(x^2)}{U(x^2)}+x\dfrac{W(x^2)}{U(x^2)}$
$U,V,W$ は、元の $P,Q$ の次数より大きくなることはない。
$x^N$ の項は、$N$ が偶数なら $\dfrac{V(x^2)}{U(x^2)}$ 側、
奇数なら $\dfrac{W(x^2)}{U(x^2)}$ 側に存在するので、どちらか片方のみ解けばよい。
求めたいのは $(x^2)^{N/2}$ の係数となるので、実質的に $N$ は半分になる。
これを再帰的に $O(\log{N})$ 回繰り返すと 求めたいものは $\dfrac{P'(x)}{Q'(x)}$ の $x^0$ の係数となり、$P'_0/Q'_0$ として求められる。
