以下の数列があるとする。
この時、以下のことが出来る。
この2つを続けて使用することで、最初の何項かがわかっているだけの数列の、かなり大きな $N$ に対する第 $N$ 項を求められる。
イメージとしては、ユークリッドの互除法を、数列（形式的べき級数）に拡張して適用したもの。
数列の最初の何項か $A=(a_1,a_2,...,a_M)$
$A$ の特性方程式。
より正確には、長さ $L+1$ の数列 $C$ であって、全ての $L \le i \le M$ に対して以下を満たすもの。
ただし、$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$ の更新の有無が変わってるんだと思う（イマイチよくわかってない）。
Berlekamp-Massey の出力は特性方程式だったが、Bostan-Mori への入力は母関数 $\dfrac{P(x)}{Q(x)}$ の $P,Q$（の係数配列）でなくてはならない。
たとえば
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
...
説明中の数列の「○○項目」は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 を畳み込む操作）
すると、$Q(x)Q(-x)$ は、$x^1,x^3,x^5,...$ の項は全て打ち消され、$x^0,x^2,x^4,...$ についての項のみ残る。
分母が $x^{偶数}$ の項しかないので、
よって、以下のようにすると
このように変形できる。
$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$ として求められる。