はじめに

　この記事は【後編】となっており、掲題の手法をC++で実装しています。

　【前編】の記事ではその手法について詳しい説明をしています。先にご覧ください。

（これらの記事の内容は「拡張フィボナッチ数列の任意の項までの和を高速に求める方法の考案」¹を加筆、修正したものです。）

プログラムコード（C++）

　$k$-ボナッチ数列の項を求めるの計算は行列の考え方を用いています。

　フィボナッチ数列に黄金比による一般項があるのと同様に、$k$-ボナッチ数列の一般項も知られているようです²。しかし一般項には無理数が含まれ、計算する値が大きくなると誤差が生じます。そのため今回は一般項を使わず整数のみで計算しています。同じ理由で行列の対角化も行っていません。

　以下のコードが作成した $k$-ボナッチ数列に関する計算を行うクラスのコードです（個人色の濃いコードになっています...）。

Fibonacci.hpp

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

/*k-ボナッチ数列
（数列の定義は,{F(0)=F(1)=...=F(k-2)=0, F(k-1)=1, F(n)=F(n-1)+F(n-2)+...+F(n-k)}）*/
class Fibonacci{
    using vl  = vector<long long>;
    using vvl = vector<vl >;
    const int mod=1e9+7;
    const int k;        //k:=(直前の足す項の数). （2以上にする! 100以下程度が実用的?）
    vvl e,a,ia,na,f;    //e:=(k次の単位元), a:=(k次の正方行列A), ia=(Aの逆行列), na:=(A^k), f:=(k*1行列).
    int inv;            //k-1の逆数.

    void err(){cerr<<"err> There are somehing err."<<endl;}
    void build(){
        REP(i,k) e[i][i]=a[0][i]=1LL;
        REP(i,k-1) a[i+1][i]=ia[i][i+1]=1LL;
        ia[k-1][0]=1LL;
        for(int i=1;i<k;++i) ia[k-1][i]=-1LL;
        na=matrix_pow(a,k);
        f[0][0]=1LL;
        inv=mod_pow(k-1,mod-2);
    }
    int mod_pow(long long x,unsigned int n){//x^n(mod 1e9+7). O(log N).
        if(n==0) return 1;
        long long res=mod_pow(x*x%mod,n/2);
        if(n&1) res=res*x%mod;
        return (int)res;
    }
    vvl matrix_mul(vvl &x,vvl &y){//行列の積をとる. O(K^3).
        if(x[0].size()!=y.size()){
            err();
            return e;
        }
        vvl res(x.size(),vl(y[0].size(),0LL));
        REP(i,x.size())REP(j,y[0].size())REP(l,x[0].size()){
            res[i][j]=(res[i][j]+x[i][l]*y[l][j]%mod)%mod;
        }
        return res;
    }
    vvl matrix_pow(vvl x,unsigned long long n,vvl res=vvl()){//行列xの累乗と行列resの積をとる. O((K^3)*log N).
        if(res.size()==0) res=e;
        while(n>0){
            if(n&1LL) res=matrix_mul(x,res);
            x=matrix_mul(x,x);
            n>>=1;
        }
        return res;
    }

public:
    /*constructor.*/
    Fibonacci(unsigned int k_=2)
    :k(k_),e(k_,vl(k_,0LL)),a(k_,vl(k_,0LL)),ia(k_,vl(k_,0LL)),f(k_,vl(1,0LL)){
        if(k_<2) err();
        build();
    }
    template<typename Type>
    Fibonacci(vector<Type> v_)
    :k(v_.size()),e(v_.size(),vl(v_.size(),0LL)),a(v_.size(),vl(v_.size(),0LL)),
    ia(v_.size(),vl(v_.size(),0LL)),f(v_.size(),vl(1,0LL)){
        if(v_.size()<2) err();
        build();
        REP(i,k) f[k-1-i][0]=(v_[i]+mod)%mod;
    }

    int fibonacci(long long n){//k-ボナッチ数列のn項目を返す. O((K^2)*log N).
        if(0<=n and n<=k-1) return f[k-1-n][0];
        vvl ff=(n<0?matrix_pow(ia,-n,f):matrix_pow(a,n-(k-1),f));
        return (n<0?(int)ff[k-1][0]:(int)ff[0][0]);
    }
    /*k-ボナッチ数列のn項目から連続した長さlenの部分列を返す. O((K^2)*log N+(Len/K)*K^2).*/
    vector<int> fibonacci(long long n,unsigned int len){
        vector<int> res(len);
        vvl ff=(n<0?matrix_pow(ia,-n,f):matrix_pow(a,n,f));
        for(long long i=n;i<n+len;i+=k){
            for(int j=0;j<k and i+j<n+len;++j) res[i+j-n]=(int)ff[k-1-j][0];
            ff=matrix_mul(na,ff);
        }
        return res;
    }
    int sum(long long n){//k-ボナッチ数列の0項目からn項目までの和を返す. O((K^2)*log N).
        if(0<=n and n<=k-2) return 0;
        if(n==k-1) return 1;
        vector<int> mem,memm;
        if(n<0) mem=fibonacci(1LL,k), memm=fibonacci(n,k);
        else mem=fibonacci(n+1,k), memm=fibonacci(0LL,k);
        long long res=((long long)mem[k-1]-memm[k-1]+mod)%mod;
        for(int j=0;j<=k-3;++j){
            long long tmp=(j+1)*(((long long)mem[k-3-j]-memm[k-3-j]+mod)%mod)%mod;
            res=(res-tmp+mod)%mod;
        }
        res=res*inv%mod;
        return (int)res;
    }
};

概要説明

　説明のため、以降より次のようになる $k$ 次の正方行列と $k×1$ 行列をそれぞれ $A, x$ とします。

\begin{align}
A^{n}x
=\begin{bmatrix}
1 & 1 & 1 & \cdots & 1 & 1\\
1 & 0 & 0 & \cdots & 0 & 0\\
0 & 1 & 0 & \cdots & 0 & 0\\
0 & 0 & 1 & \cdots & 0 & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\
0 & 0 & 0 & \cdots & 1 & 0
\end{bmatrix}^{n}
\begin{bmatrix}
F_{k-1}\\
F_{k-2}\\
\vdots\\
F_2\\
F_1\\
F_0
\end{bmatrix}
=\begin{bmatrix}
F_{n+k-1}\\
F_{n+k-2}\\
\vdots\\
F_{n+2}\\
F_{n+1}\\
F_{n}
\end{bmatrix}\\
\end{align}

　$e, a, ia, na, f$ は行列用の2次元配列です。$e$ は $k$ 次の単位行列、$a$ は行列 $A$、$ia$ は $A^{-1}$、$na$ は $A^k$、$f$ は行列 $x$ を表します。

　$inv$ は $k-1$ の逆元です³。

　$k$-ボナッチ数列は項が大きくなるにつれ数の値が非常に大きくなるので、素数である $10^9+7$ で $mod$ をとっています。

コンストラクタ

$Fibonacci(\ k\ )$：$k$-ボナッチ数列に関した計算をするオブジェクトを定義する.
$Fibonacci(\ v\ )$：配列 $v$ に格納された数をはじめの項に設定した数列を定義する.

　オブジェクトを定義するときの引数は$k$ の値、または$0$ 項目から $k-1$ 項目の値を格納した配列です。$vector<int>(\{2,1\})$ を引数として渡すとリュカ数列⁴がつくれます。

　$k$ の値が2未満だとバグるります（エラー）。

その他

$err()$：エラー警告を行う.
$build()$：計算に使う行列や整数の数値を初期化する.
$mod\_pow( x, n )$：$x^n(mod\ 10^9+7)$を計算する. 繰り返し二乗法を利用. 計算量は$O(\log n)$.
$matrix\_mul( x, y )$：行列 $x$ と $y$ の積を計算する. 最悪計算量 $O(k^3)$.
$matrix\_pow( x, n, res )$；行列 $x^n$ と $res$ の積を計算する. 繰り返し二乗法を利用. 最悪計算量 $O(k^3 \log n)$.

　行列の積の計算では、例えば$a×b$ 行列と $b×c$ 行列の積をとる場合、$a×b×c$ 回計算を行います。そのため関数 $matrix\_pow( x, n, res )$ では、繰り返し二乗法の途中で第3引数の行列に直接かけさせることにより（必要な $c$ の値を最小限にする）、可能な限り計算回数を減らす工夫をしました。

機能の説明

　利用できる主な機能は「任意の項の値を求める」「任意の項までの和を求める」の2つです。

数列の任意の項を求める

$fibonacci( n )$：$k$-ボナッチ数列の $n$ 項目の数を返す. $O(k^3 \log n)$.
$fibonacci( n, len )$：$k$-ボナッチ数列の $n$ 項目から長さlenの部分列を返す. $O((k^3 \log n)+\frac{len}{k}×k^2)$.

　実は負の項への計算にも拡張してます。このとき $A$ の逆行列を利用します。

\begin{align}
(A^{-1})^{n}x
=\begin{bmatrix}
0 & 1 & 0 & 0 & \cdots & 0\\
0 & 0 & 1 & 0 & \cdots & 0\\
0 & 0 & 0 & 1 & \cdots & 0\\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\
0 & 0 & 0 & 0 & \cdots & 1\\
1 & -1 & -1 & -1 & \cdots & -1
\end{bmatrix}^{n}
\begin{bmatrix}
F_{k-1}\\
F_{k-2}\\
\vdots\\
F_2\\
F_1\\
F_0
\end{bmatrix}
=\begin{bmatrix}
F_{-n+k-1}\\
F_{-n+k-2}\\
\vdots\\
F_{-n+2}\\
F_{-n+1}\\
F_{-n}
\end{bmatrix}\\
\end{align}

　また、行列$A^k$をあらかじめ計算することにより、1回の行列の積の計算で次の $k$ 個の項の値を求められるように工夫しています。

\begin{align}
A^{k}x
\begin{bmatrix}
1 & 1 & \cdots & 1 & 1\\
1 & 0 & \cdots & 0 & 0\\
0 & 1 & \cdots & 0 & 0\\
\vdots & \vdots & \ddots & \vdots & \vdots\\
0 & 0 & \cdots & 1 & 0
\end{bmatrix}^{k}
\begin{bmatrix}
F_{k-1}\\
F_{k-2}\\
\vdots\\
F_1\\
F_0
\end{bmatrix}
=\begin{bmatrix}
F_{2k-1}\\
F_{2k-2}\\
\vdots\\
F_{k+1}\\
F_{k}
\end{bmatrix}\\
\end{align}

数列の和を求める

$sum( n )$：数列の0項目からn項目までの数の和を求める. 計算量は$O(k^3 \log n)$.

　次の方程式を用いています。式の導出は【前編】をご参照ください。

\begin{align}
\sum^n_{i=0}F_i=\frac{F_{n+k}-F_{k-1}- \begin{align} \sum^{k-3}_{j=-1}\Bigl\{ (j+1)(F_{n+k-2-j}-F_{k-3-j})\Bigr\} \end{align} }{k-1}\\\\
\text(ただし, k\geq2かつn\geq0)
\end{align}

\begin{align}
\sum^n_{i=0}F_i=\frac{F_{k}-F_{n+k-1}- \begin{align} \sum^{k-3}_{j=-1}\Bigl\{ (j+1)(F_{k-2-j}-F_{n+k-3-j})\Bigr\} \end{align} }{k-1}\\\\
\text(ただし, k\geq2かつn<0)
\end{align}

　必要な項は長さ $k$ の部分数列が2つのみです。ですので、従来法よりかなり高速に求められることがわかります。

あとがき

　前編、後編あわせて以上で終わりです。

　フィボナッチ数列も調べてみるとなかなか不思議ですね。

（未熟者ですので、記事の内容に誤りなどがあるかもしれません。もし何か気付かれた方がいらっしゃいましたら、コメントなどで知らせていただけると有難いです。）

today Diary「拡張フィボナッチ数列の任意の項までの和を高速に求める方法の考案」（はてなブログ）. ↩
@magolorsさんの記事「拡張フィボナッチ数列(k-ボナッチ数列 : 直前のk項の和)の一般項を線形代数とPythonで求める」を参照. ↩
$x(mod\ p)$ の逆元は $x^{p-2}(mod\ p)$ です. 詳しくは@drkenさんの記事「フェルマーの小定理の証明と使い方」を参照. ↩
Wikipwdiap「リュカ数」を参照. ↩

拡張フィボナッチ数列（k-ボナッチ数列）の任意の項までの和を高速に求める手法の考案【後編】