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

とりあえず使うことを目的としたメモで、原理の解説はしていない。

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

  • 無限に続き、最初の何項か $a_1,a_2,\ldots,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$ 項を求められる。

Barlekamp-Massey algorithm

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

入力

数列の最初の何項か $A=(a_1,a_2,\ldots,a_M)$

出力

$A$ の特性方程式。

より正確には、長さ $L+1$ の数列 $C$ であって、全ての $L \le i \le M$ に対して以下を満たすもの。

  • $C_0a_i+C_1a_{i-1}+C_2a_{i-2}+\ldots+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$ の更新の有無が変わってるんだと思う(イマイチよくわかってない)。

Bostan-Mori algorithm

入力

$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})$ となる。

programming_algorithm/number_theory/barlekamp_massey.txt · 最終更新: 2021/04/17 by ikatakos
CC Attribution 4.0 International
Driven by DokuWiki Recent changes RSS feed Valid CSS Valid XHTML 1.0