数列の漸化式の特定と第N項の計算

以下の数列があるとする。

• 無限に続き、最初の何項か $a_1,a_2,...,a_M$ がわかっている
• あんまり多くない次数（具体的にわからなくてもいい）の線形漸化式で表現できる
  • ○ $a_n = 3a_{n-1} + 5a_{n-2} + 7a_{n-3}$
  • × $a_n = a_{n-1}^2$
  • × $a_n = n \times a_{n-1}$

この時、以下のことが出来る。

• Barlekamp-Massey（バーレカンプ-マッシー）のアルゴリズム
  • 特性方程式を特定する
    • 漸化式を $a_n = 3a_{n-1} + 5a_{n-2}$ とすると、$x^2 - 3x - 5 = 0$ みたいになるやつ
    • ここから漸化式の形もすぐ求められる
    • あり得るもののうち、次数が最小のものが求まる
      • そのため、はじめに与える数列の長さが足りないと、本来の特性方程式が求まる保証は無い
      • 正解が $d$ 次漸化式の場合、$M=2d$ 項あれば十分
• Bostan-Mori のアルゴリズム
  • 線形漸化式で表される数列の第 $N$ 項を高速に計算する
  • $d$ 次漸化式として、$O(d^2 \log{N})$ または $O(d \log{d} \log{N})$

この2つを続けて使用することで、最初の何項かがわかっているだけの数列の、かなり大きな $N$ に対する第 $N$ 項を求められる。

イメージとしては、ユークリッドの互除法を、数列（形式的べき級数）に拡張して適用したもの。

• Berlekamp–Massey_algorithm
• 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$

英語版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$ の更新の有無が変わってるんだと思う（イマイチよくわかってない）。

• 線形漸化式を満たす数列の N 項目を計算するアルゴリズム

$A$ を、$d$ 次線形漸化式で表される数列とする。

①「正整数 $N$」と、

②「$A$ の母関数が $G(x)=\dfrac{P(x)}{Q(x)}$ と表される時の、$d$ 次式 $P$ の係数および $d+1$ 次式 $Q$ の係数」

が入力となる。 ただし②は③からも計算できる。

③「$A$ の漸化式の各係数、および最初の $d$ 項」

たとえば
a[1] = a[2] = a[3] = 1
a[n] = a[n-1] + 2*a[n-2] + 3*a[n-3]  のとき

③
最初の3項: [1, 1, 1]
Aの係数  : [1, 2, 3]
↓
②
Aの特性方程式の係数をQとする（ただし第1項は 1）
  Q: [1, -1, -2, -3]
Qと最初の3項を（たたみ込みの容量で）かけ、最初の3項をPとする
  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
                                  ...

Barlekamp-Masseyからの引き継ぎでいうなら、出力がそのまま $Q$ となっていて、 $A$ の最初の何項かはBarlekamp-Masseyの入力でもあるので、かけあわせて $P$ が求められる。

数列 $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] が答え

たたみ込みを愚直にやると $O(d^2 \log{N})$、高速フーリエ変換などを使えるなら $O(d\log{d}\log{N})$ となる。