拡張フィボナッチ数列の任意の項までの和を高速に求める方法の考案

 暑い日が続いていますね...。最近2025大阪万博のロゴが気になってます。


 今回は掲題の内容について、高校数学の知識から実用的な(?)プログラムができたので紹介したいと思います。

 少し長くなりますがお付き合いください。

(未熟者ですので細かい点において誤った情報になっている可能性があります。もしお気づきの方がいらっしゃいましたら、ご指摘していただけると助かります。)

まえがき

 フィボナッチ数列とは「前にある2つの数を足した数の列」で、

 
\begin{cases}
f_0 = 0,  f_1 = 1\\
f_k = f_{k-1} + f_{k-2}  (k \geq 2)
\end{cases}

と定義されます。

 他にも前3つの数を足していくトリボナッチ数列や、前4つを足していくテトラボナッチ数列などがあります。


 ところでフィボナッチ数列 k番目の項までの和については、


 
\begin{eqnarray}
\sum^k_{i=0} f_i = f_{k+2} - f_1 = f_{k+2} - 1
\end{eqnarray}


という1つの一般項から導出できることが知られています。この式の求め方についてはこちらの他の方のブログをご覧ください。


 私は先日、
「この式の考え方を応用したらトリボナッチ数列やテトラボナッチ数列でも同じような式が書けるんじゃない?!」
という発想に至り、一般化した式を導出しました。Σ(゚Д゚)スゲェ!!



概要 ~数列の定義と“和の式”~

 以下の文章から足す数の個数を n, その数列を A(n)=\{a_0, a_1, a_2, ......\}とします。

拡張フィボナッチ数列 A(n)の定義

 すると数列 A(n)は、

 
\begin{cases}
a_0 = a_1 = ...... = a_{n-2} = 0\\
a_{n-1} = 1\\
a_k = a_{k-1} + a_{k-2} + ...... + a_{k-n}   (k \geq 2)
\end{cases}

と定義できます。

数列 A(n)の任意の項までの和の式

 このとき、数列 A(n) k番目の項までの和は



\begin{eqnarray}

\sum^k_{i=0} a_i = \frac{ a_{k+n} - a_{n-1} -\begin{eqnarray} \sum^{n-3}_{j=0} \{ (j+1)(a_{k+n-2-j} - a_{n-3-j}) \}\end{eqnarray} }{n-1}

\end{eqnarray}


と表現できます。これが今回私が導出した式です!

 この式のすごいところは、 k+1個の一般項の和の式が、僅か 2+2(n-2)個の一般項のみから導き出せる」ということです!



式の導出と証明

 タイプが面倒くさかったので自分が導出したときのメモを張っておきます(自我乱れて見づらいですが...)

 数学的帰納法で証明しました。(途中 n=3, 4, 5について書いていますが蛇足なので読み飛ばしてください。)




完成したプログラムFibonacci(C++

 以下が作ったプログラムです(個人色強いプログラムですがご了承ください)。オーダー記法も自信がないので参考程度にお願いします。

#include <vector>
#define REP(i,n) for(int i=0;i<(n);++i)

/*フィボナッチ数列
(数列の定義は,{F(0)=...=F(n-2)=0, F(n-1)=1. F(k)=F(k-n)+F(k-(n-1))+...+F(k-2)+F(k-1).})*/
class Fibonacci{
	using llong = long long;
	using vl    = vector<llong>;
	using vvl   = vector<vl >;
	const int mod=1e9+7;
	const int n;	//n:=(前にある足す数の個数). (2以上にすること! 約100以下が実用的?)
	vvl a,b;		//a:=(n*n行列A), b:=(1*n転置行列(F(n-1),.....,F(0))').
	
	vvl mul(vvl &x,vvl &y){//n*n行列xとyの積をとる. O(n^3).
		vvl res(x.size(),vl(y[0].size(),0LL));
		if(x[0].size()==y.size()){
			REP(i,x.size())REP(j,y[0].size())REP(k,x[0].size()){
				res[i][j]=(res[i][j]+x[i][k]*y[k][j]%mod+mod)%mod;
			}
		}
		return res;
	}
	llong mod_pow(llong n,llong k){//n^k(mod p). O(log K).
		if(k==0LL) return 1LL;
		llong res=mod_pow(n*n%mod,k/2);
		if(k&1LL) res=res*n%mod;
		return res;
	}
	
public:
	Fibonacci(unsigned int n_=2):n(n_),a(n,vl(n,0LL)),b(n,vl(1,0LL)){//constructor.
		REP(i,n) a[0][i]=1LL;
		REP(i,n-1) a[i+1][i]=1LL;
		b[0][0]=1LL;
	}
	
	int fibonacci(unsigned long long k){//フェボナッチ数列のk項目を求める. O((n^3)*log(i)).
		if(k<n-1) return 0;
		if(k==n-1) return 1;
		vvl aa=a, res=b;
		while(k>0LL){
			if(k&1LL) res=mul(aa,res);
			aa=mul(aa,aa);
			k>>=1;
		}
		return res[n-1][0];
	}
	/*フェボナッチ数列のk項目までの和を求める. O((n^4)*log(i)).
	(Σ(0~k)F(i)=[F(k+n)-F(n-1)-Σ(0~n-3){(j+1)(F(k+n-2-j)-F(n-3-j))}]/(n-1)より.)*/
	int sum(unsigned long long k){
		llong res=((llong)fibonacci(k+n)-fibonacci(n-1)+mod)%mod;
		llong tmp=0LL;
		for(int j=0;j<=n-3;++j){
			llong tmpp=(j+1LL)*(((llong)fibonacci(k+n-2-j)-fibonacci(n-3-j)+mod)%mod)%mod;
			tmp=(tmp+tmpp)%mod;
		}
		res=(res-tmp+mod)%mod;
		res=res*mod_pow(n-1,mod-2)%mod;//逆元.
		return (int)res;
	}
};

構造体Fibonacciの説明

前説明

 基本的なアイデアは行列を用いています。

 項が増えると数がとても大きくなるので、素数の1,000,000,007でmodとっています。

  nの値の定義はコンストラクタの部分で行います(2未満だとバグります)。このプログラムではサイズ n^2の配列が必要となったり、 O(n^3)の計算があったりするので大きすぎる値は注意が必要です。

機能説明

 この構造体の主な機能は、数列 A(n)における

  1. 任意の項の数の取得( 関数fibonacci() ).
  2. 任意の項までの和の取得( 関数sum() ).

の2つです。
 

 関数fibonacci()は行列の考え方を利用しています(Wikipedia)。行列の積の計算では最悪 O(n^3)となりますが、繰り返し二乗法のような(?)考え方を使うので全体的に O(n^3*log (k))です(要再考!)。

 関数sum()は今回考案した考え方を利用して計算しています。sum()内ではfibonacci()を 2+2(n-2)回呼出しです。ですので計算量は大体 O(n^4*log (k))となります(要再考!)。

 つまり、 nが大きかったり項が小さい場合は愚直に調べたほうがいいかもしれません。が、求めたい項が大きい場合は有効です。



メリット・有用性

 今回考案した式の有効な場面について再度説明します。

 もし 0項目から k 項目まで数を求めて足していくということをすると、計算量は O(k * n^3 * log(k))となります(要再考!)。

 対して今回の計算量は O(n^4*log (k))です(要再考!)。

 つまり、「 nが小さく、 kが大きい」というときにより早く計算できるということになります。

実行時間の検証

 次の2つの方法において、任意の項までの和を求めるのにかかる実行時間を計測して比較します。

  1.  0項目から k 項目まで数を求めて足していく.
  2. 考案した式を用いる.

 私が試した結果は以下のようになりました。(単位は[ms]です。)

 n  k(項番号) 手法1(愚直) 手法2(新)
2 100 30.0 0.0
3 100 50.3 0.0
10 100 319.2 50.3
100 1000 (not) 354953.2
2 10000 2334.9 0.0
3 10000 46009.3 0.0
10 10000 (not) 0.0

 初めて比較計測したので正しくできているかはわかりませんが、実感でも新手法が早かったと思います。



参考文献

あとがき

 正直これ正しいのかという不安や実用性の懸念もありますが、久しぶりにイチから(いろいろ参考にしながらも)プログラムを作るのが楽しかったです。

 また冒頭でも述べたように、誤った記述があるかもしれません(もしかしたらこの記事全体が間違っているかもしれません)。もし何か気付いたという方がいらっしゃいましたら、コメント機能やツイッターにて知らせていただけると有難いです。

 あと記事を引用などする場合はリンクやタイトルなどの明示をお願いいたします。