拡張されたフィボナッチ数列をPythonで実装してみるPart2 (項数の変更)

はじめに

フィボナッチ数列をPythonで実装するシリーズ記事のPart2です。Part0とPart1を読了したことを前提としているのでまだの方は是非お読みください。

• Part0:Pythonによるフィボナッチ数列の色々な求め方
• Part1:負の数番を考える
• Part2:項数の一般化 (イマココ)
• Part3:初期値の一般化

ではPart2の今回はフィボナッチ数列の項数をPython3で変更していきましょう。

項数の変更

フィボナッチ数列は各項が前二項の足し算で表される数列です。
ならば前三項の足し算で表される数列や前四項の足し算で表される数列も同様に考えることができそうです。
このような項数の変更について考えていきましょう。

トリボナッチ数列、テトラナッチ数列

前三項、前四項の足し算で表される数列をそれぞれトリボナッチ数列、テトラナッチ数列といいます。それぞれの定義を確認してみましょう。

トリボナッチ数列は非負の整数nについて三つの初期条件(初項$T_0=0$,第二項$T_1=0$,第三項$T_2=1$)を持つ漸化式$T_n=T_{n-1}+T_{n-2}+T_{n-3} $で表される数列です。具体的には$0, 0, 1, 1, 2, 4, 7, 13, 24, 44, 81, 149, 274,\cdots$ です。

テトラナッチ数列は非負の整数nについて四つの初期条件(初項$T_0=0$,第二項$T_1=0$,第三項$T_2=0$,第四項$T_3=1$)を持つ漸化式$T_n=T_{n-1}+T_{n-2}+T_{n-3}+T_{n-4} $で表される数列です。具体的には$0, 0, 0, 1, 1, 2, 4, 8, 15, 29, 56, 108, 208, 401,\cdots$ です。

これら二つの数列を Part0 で示した3つの方法で$T_n$を求める関数Trib(n)とTetra(n)をそれぞれ実装してみます。

反復法

変数を増やすだけです。

trib_tetra_repetition.py

def Trib(n):
    a, b, c = 0, 0, 1
    if n == 0:
        return a
    elif n == 1:
        return b
    elif n == 2:
        return c
    else:
        for i in range(n-2):
            a, b, c = b, c, a + b + c
        return c

print([Trib(n) for n in range(15)])
# [0, 0, 1, 1, 2, 4, 7, 13, 24, 44, 81, 149, 274, 504, 927]

def Tetra(n):
    a, b, c, d = 0, 0, 0, 1
    if n == 0:
        return a
    elif n == 1:
        return b
    elif n == 2:
        return c
    elif n == 3:
        return d
    else:
        for i in range(n-3):
            a, b, c, d = b, c, d, a + b + c + d
        return d

print([Tetra(n) for n in range(15)])
# [0, 0, 0, 1, 1, 2, 4, 8, 15, 29, 56, 108, 208, 401, 773]

出力は成功しているようです。ここではフィボナッチ数列のときの実装をそのまま拡張しましたが、辞書型を使えばもっとすっきり書けますね。

再帰法

項数を増やすだけです。

trib_tetra_recursion.py

def Trib(n):  
    if n in [0,1]:
        return 0
    elif n == 2:
        return 1
    else:
        return Trib(n-1) + Trib(n-2) + Trib(n-3)

print([Trib(n) for n in range(15)])
# [0, 0, 1, 1, 2, 4, 7, 13, 24, 44, 81, 149, 274, 504, 927]

def Tetra(n):  
    if n in [0,1,2]:
        return 0
    elif n == 3:
        return 1
    else:
        return Tetra(n-1) + Tetra(n-2) + Tetra(n-3) + Tetra(n-4)
    
print([Tetra(n) for n in range(15)])
# [0, 0, 0, 1, 1, 2, 4, 8, 15, 29, 56, 108, 208, 401, 773]

出力は成功しているようです。漸化式をそのまま書くだけなのでこれがやはり最も可読性が高いですね。

一般項

Part0の一般項の求解 で示した通り、sympyには漸化式を解く機能があります。これを使えばトリボナッチ数列やテトラナッチ数列の一般項を求めることができそうです。
やってみましょう。

.py

import sympy as sym
f = sym.simplify("f(n)") 
s = sym.simplify("f(n)-f(n-1)-f(n-2)-f(n-3)")
ini = sym.simplify({0:0,1:0,2:1})

print(sym.rsolve(s,f,ini)) #初期条件iniの漸化式sをfについて解く
# NotImplementedError: multiple candidates for the minimal polynomial of -38*(3*sqrt(33) + 19)**(1/3) - 6*sqrt(33)*(3*sqrt(33) + 19)**(1/3) - 8*(3*sqrt(33) + 19)**(2/3) - 6*sqrt(33)

失敗しました。エラーコードを読む限り、「解の最小多項式に複数の候補があるからsympyに実装されている機能では解けないよ」ということらしいです。トリボナッチ数列の漸化式は直接解くには複雑すぎるということでしょうか。
テトラナッチ数列についても同様に試みましたが、私の環境ではフリーズしてしまい、エラーメッセージすら見ることができませんでした。

しかし、この程度では諦めません。アプローチの方法を変えましょう。
Part0の補足 で、行列を使った解き方を紹介しました。あれを使います。
前回記事では湧いて出た定理のように扱いましたが、実はあの行列の関係式はフィボナッチ数列の一般項を導出する際に出てくる関係式なのです。説明していきます。
フィボナッチ数列の漸化式$F_n=F_{n-1}+F_{n-2}$は行列の積で

F_n=\left(\begin{array}{rrr}1 & 1 \end{array}\right)\left(\begin{array}{rrr} F_{n-1}  \\F_{n-2} \end{array}\right)

と表記できます。
また、当たり前ですが

F_{n-1}=\left(\begin{array}{rrr}1 & 0 \end{array}\right)\left(\begin{array}{rrr} F_{n-1}  \\F_{n-2} \end{array}\right)

です。
これらを使って帰納的に以下のことがいえます。

\left(\begin{array}{ccc} F_{n}  \\F_{n-1} \end{array}\right) =
 \left(\begin{array}{rrr}1 & 1  \\1 & 0\end{array}\right)
\left(\begin{array}{rrr} F_{n-1}  \\F_{n-2} \end{array}\right)=
 \left(\begin{array}{rrr}1 & 1  \\1 & 0\end{array}\right)^2
\left(\begin{array}{rrr} F_{n-2}  \\F_{n-3} \end{array}\right)=
\cdots
=
 \left(\begin{array}{rrr}1 & 1  \\1 & 0\end{array}\right)^{n-1}
\left(\begin{array}{rrr} F_{1}  \\F_{0} \end{array}\right)=
 \left(\begin{array}{rrr}1 & 1  \\1 & 0\end{array}\right)^{n-1}
\left(\begin{array}{rrr} 1  \\ 0 \end{array}\right)

ここで~\left(\begin{array}{rrr}1 & 1  \\1 & 0\end{array}\right)^{n-1}=
\left(\begin{array}{rrr}A & B  \\C & D\end{array}\right)~とすると
\left(\begin{array}{ccc} F_{n}  \\F_{n-1} \end{array}\right)=
\left(\begin{array}{rrr}A & B  \\C & D\end{array}\right)\left(\begin{array}{rrr} 1  \\ 0 \end{array}\right)~であるから
F_n=A~が成り立つ。

これでPart0で示した定理が証明できました。この行列の関係式からジョルダン標準形などの数学的な手続きを踏んで一般項を導出するのですが、sympyでは行列の計算が簡単にできるので、この行列の関係式が一般項であるといっても差支えないでしょう。

同様の議論でトリボナッチ数列、テトラナッチ数列についても以下のことが言えます。
・トリボナッチ数列

\left(\begin{array}{ccc} T_{n}  \\T_{n-1} \\T_{n-2} \end{array}\right) =
 \left(\begin{array}{rrr}1 & 1 & 1 \\1 & 0 & 0\\ 0 & 1 & 0\end{array}\right)
\left(\begin{array}{ccc} T_{n-1}  \\T_{n-2} \\T_{n-3} \end{array}\right)=
\cdots
=
 \left(\begin{array}{rrr}1 & 1 & 1 \\1 & 0 & 0\\ 0 & 1 & 0\end{array}\right)^{n-2}
\left(\begin{array}{ccc} T_{2}  \\T_{1} \\T_{0} \end{array}\right)=
 \left(\begin{array}{rrr}1 & 1 & 1 \\1 & 0 & 0\\ 0 & 1 & 0\end{array}\right)^{n-2}
\left(\begin{array}{ccc} 1  \\0 \\0 \end{array}\right)

ここで~ \left(\begin{array}{rrr}1 & 1 & 1 \\1 & 0 & 0\\ 0 & 1 & 0\end{array}\right)^{n-2}=
 \left(\begin{array}{rrr}A & B & C \\D & E & F\\ G & H & I\end{array}\right)~とすると
T_n=A~が成り立つ。

・テトラナッチ数列

 \left(\begin{array}{rrr}1 & 1 & 1 & 1 \\ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{array}\right)^{n-3}=
 \left(\begin{array}{rrr}A & B & C & D \\ E & F & G & H \\ I & J& K & L \\ M & N & O & P \end{array}\right)~とすると
T_n=A~が成り立つ。

この一般項(行列の関係式)を使って$T_n$を求める関数Trib(n)とTetra(n)をそれぞれ実装してみましょう。

.py

def Trib(n):
    if n == 0 or n == 1:
        return 0
    else:
        A = sym.Matrix([
            [1,1,1],
            [1,0,0],
            [0,1,0]
        ])
        FibArray=A**(n-2)
        return FibArray[0,0]

print([Trib(n) for n in range(15)])
# [0, 0, 1, 1, 2, 4, 7, 13, 24, 44, 81, 149, 274, 504, 927]
   
def Tetra(n):
    if n == 0 or n == 1 or n == 2:
        return 0
    else:
        A = sym.Matrix([
            [1,1,1,1],
            [1,0,0,0],
            [0,1,0,0],
            [0,0,1,0]
        ])
        FibArray=A**(n-3)
        return FibArray[0,0]
    
print([Tetra(n) for n in range(15)])
# [0, 0, 0, 1, 1, 2, 4, 8, 15, 29, 56, 108, 208, 401, 773]

どちらも成功です。数学的な仕組みさえわかれば簡単ですね。

ここでは一般項を導くために数学的なプロセスをいくつか経ましたが、トリボナッチ数列orテトラナッチ数列を一般項を用いて求めたいだけならWikipedia等に一般項の式が載っているのでそれを使うのが簡単です。しかし、これ以降の拡張のことを考えて、あえて行列による解法を選びました。

項数に応じた名称について

本題は少し休憩して、余談です¹。

フィボナッチ数列、トリボナッチ数列、テトラナッチ数列、これらの名称について、Wikipediaを読むとこんなことが書いてありました。

これらフィボナッチ数列の類似物を、項数 k に対応するラテン語またはギリシャ語に由来する倍数接頭辞を「フィボナッチ」と組み合わせた名称で呼ぶ^{[注釈 1]}。

注釈1:当然のことだが "Fibonacci" は人名であって、"fibo-" + "-nacci" や "fi-" + "-bonacci" という構成の合成語でもないし、もちろん "fi-" や "fibo-" が "2" の意味を持つわけでもない（ただし、摩擦音 f と破裂音 b が音韻的に近い関係にあることから 2 を表す "bi-" を "fi-" に結び付けての類推ではあるかもしれない）が、「フィボナッチ」の語を頭から適当な音節分だけ倍数を表す接頭辞で置き換えるという、冗談のような名付けになっている。

こういう数学者たちの遊び心みたいなの、いいですよね。

さて、次にお話しするのは項数を一般化したフィボナッチ数列の類似数列です。ここではトリボナッチやテトラナッチに倣い、n-ボナッチ(エヌボナッチ)数列と呼ぶことにしましょう。

n-ボナッチ数列

n-ボナッチ数列を非負の整数$n,k ~(n\geqq2)$ についてn個の初期条件($N_0=N_1=\cdots=N_{n-2}=0$ , $N_{n-1}=1$)を持つ漸化式$N_k=N_{k-1}+N_{k-2}+\cdots+N_{k-n} $で表される数列であると定義します²。
つまり、フィボナッチ数列は$n=2$のときのn-ボナッチ数列であるといえます。

ではn-ボナッチ数列を実装する方法を考えていきましょう。

実装ルール

実装する前に、実装ルールを決めておきます。
これから考えるn-ボナッチ数列の$N_k$を求める自作関数は引数を二つ持ちます。一つ目はn-ボナッチ数列のnにあたる引数で、もう一つは$N_k$の$k$を指定する引数です。n_bonacci(n,k)のように実装します。つまり、n_bonacci(2,8)ならば、フィボナッチ数 $F_8$ (すなわち21)を戻り値として出力します。

以下に再帰法と一般項を使う方法の二つの実装例を紹介します。反復法でももちろんできますが、実装方法は再帰法と大差ないので省きました。興味のある方は是非、反復法で実装してコメントに書き込んでください。喜びます。

再帰法

再帰法は漸化式をそのまま実装するだけなのでそこまで複雑にはなりません。ただ、項数が定数ではないので、少し工夫が必要です。for文を使いましょう。
n-ボナッチ数列の漸化式$N_k=N_{k-1}+N_{k-2}+\cdots+N_{k-n} $はfor文を用いて次のように表せます。

.py

# n,kは与えられているとする
recurrence_relation = 0 #漸化式の初期化
for i in range(1, n+1):
    recurrence_relation += n_bonacci(n, k-i)

これを使って実装していきましょう。

.py

def n_bonacci(n,k):
    if k <= n-2:
        return 0
    elif k == n-1:
        return 1
    else:
        recurrence_relation = 0
        for i in range(1,n+1):
            recurrence_relation += n_bonacci(n, k-i)
        return recurrence_relation
    
print([n_bonacci(2, k) for k in range(15)])  #フィボナッチ数列
# [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377] 

print([n_bonacci(3, k) for k in range(15)])  #トリボナッチ数列
# [0, 0, 1, 1, 2, 4, 7, 13, 24, 44, 81, 149, 274, 504, 927]

print([n_bonacci(4, k) for k in range(15)])  #テトラナッチ数列
# [0, 0, 0, 1, 1, 2, 4, 8, 15, 29, 56, 108, 208, 401, 773]

成功しているようです。n=8にしたときでもオンライン整数列大辞典の オクタナッチ数列 と一致していました。

一般項

トリボナッチ数列以降はsympyで漸化式(階差方程式)を直接求めることができないので、行列を使います。
トリボナッチ数列、テトラナッチ数列と同様の議論でn-ボナッチ数$N_k~(k\geqq n)$について以下のことが言えます。

n \times n 行列~A= \left(\begin{array}{ccc}
1 & 1 & \cdots & \cdots &\cdots & 1 \\
1 & 0 & \cdots & \cdots &\cdots & 0 \\ 
0 & 1 & \ddots &  & & \vdots \\ 
\vdots & \ddots & \ddots &\ddots& & \vdots  \\
\vdots &  & \ddots& \ddots & \ddots & \vdots  \\
0&\cdots&\cdots&0&1&0 \\
\end{array}\right)
について~A^{k-(n-1)}の一行一列要素はN_kと等しい。

まずはこの$n \times n $行列をつくるロジックを考えます。この行列は第一行は全て$1$、第$n$行$(n\geqq 2)$は$n-1$列目だけが$1$で他は全て$0$であることが分かっています。
入れ子構造のリスト内法表記で行列を全部一気に作るのも良いですが、行列の結合を使ってみましょう。
この$n \times n $行列は以下の三つの行列を結合したものです。

(n-1) \times (n-1) 単位行列~\left(\begin{array}{ccc}
1 & 0 & \cdots & \cdots & 0  \\
0 & 1 & \ddots &  &\vdots  \\ 
\vdots & \ddots & \ddots & \ddots &\vdots  \\ 
\vdots &  &\ddots &\ddots & 0 \\
0 & \cdots & \cdots& 0 & 1  \\
\end{array}\right)
~,~
1\times n 行列~\left(\begin{array}{ccc}
1 &  \cdots & 1  \\
\end{array}\right)
~,~
(n-1)\times 1 行列\left(\begin{array}{ccc}
0 \\  \vdots \\ 0  
\end{array}\right)

$n \times n $単位行列はsym.eye(n)、AとBの行結合はA.row_join(B)、AとBの列結合はA.col_join(B)と書けます。

.py

import sympy as sym
def n_nMatrix(n):
    A1 = sym.Matrix([[0] for i in range(n-1)]) #(n-1)×1行列
    A2 = sym.eye(n-1).row_join(A1) #(n-1)×(n-1)単位行列と行結合
    A = sym.Matrix([[1 for i in range(n)]]).col_join(A2) #1×n行列と列結合
    return A   
     
print(n_nMatrix(4))
# Matrix([[1, 1, 1, 1], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])

これで一般項の導出に必要な行列はつくれました。
では、n-ボナッチ数列の$N_k$を求める自作関数n_bonacci(n,k)を実装してみましょう。

.py

import sympy as sym
def n_bonacci(n,k):
    if k <= n-2:
        return 0
    else:
        A1 = sym.Matrix([[0] for i in range(n-1)])
        A2 = sym.eye(n-1).row_join(A1) 
        A = sym.Matrix([[1 for i in range(n)]]).col_join(A2) 
        n_bonacciMatrix=A**(k-(n-1))
        return n_bonacciMatrix[0,0] #1行1列要素を返す

print([n_bonacci(2, k) for k in range(15)])  #フィボナッチ数列
# [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377] 

print([n_bonacci(3, k) for k in range(15)])  #トリボナッチ数列
# [0, 0, 1, 1, 2, 4, 7, 13, 24, 44, 81, 149, 274, 504, 927]

print([n_bonacci(4, k) for k in range(15)])  #テトラナッチ数列
# [0, 0, 0, 1, 1, 2, 4, 8, 15, 29, 56, 108, 208, 401, 773]

おわりに

今回はフィボナッチ数列の項数を拡張してn-ボナッチ数列を実装しました。

Part3では初期値の変更を行います。

誤字脱字や誤った記述、更に良いアルゴリズム、気づいたこと、分かりにくい箇所、感想などがあればぜひコメントをください。

参考文献

フィボナッチ数 - Wikipedia

1. 「閑話休題」という語を使おうと思いましたが、調べてみると誤用でした。本来、閑話休題は「無駄話をやめて本題に戻ること」を意味するのだそうです(そして逆の意味に誤用されがちだそうです)。勘違いしていました。 ↩

2. この記事でこのように定義するだけで、一般的な名称、定義ではありません。悪しからず。 ↩