フィボナッチ数列で覚えるkitamasa法

最近kitamasa法というものを覚えたので解説する
フィボナッチ数列ではあまりうまみがないが具体例ということでフィボナッチ数列を使用する

# フィボナッチ数列
フィボナッチ数列は前の二つの数字を足したものが次の項になる、という数列で以下の漸化式で定義されます。

```math
F_0 = 0, \\
F_1 = 1, \\
F_{n+2} = F_n + F_{n+1} (n ≥ 0)
```
実際には以下のような
数値になります(第0～21項)

```
0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946,
```

# 愚直な計算方法
前二つの数字を覚えておきさえすればいいので項の小さいほうから計算すれば容易に計算できます
この場合計算回数は第$N$項を求めるのに$N$回の計算が必要です。

# 行列累乗
[行列累乗](https://qiita.com/flour/items/9f8a20112f45f0ee5f71)という方法でも計算ができます 
この方法だとフィボナッチ数列の場合第$N$項を求めるのに$log(N)$回ほどで済みます。
ただし、行列の積は $ n \times n $ の行列の場合 $ n^3 $ かかります。
フィボナッチ数列の場合は $ 2 \times 2 $ の行列なのでよいのですが前$K$項を足すフィボナッチ数列の一般化を考えたときに計算量が $O(K^3 logN)$ になってKが1000などの場合は計算に時間がかかるようになります。
# kitamasa法

ここでkitamasa法という方法を使うと前$K$項を足すフィボナッチの場合の計算量 $O(K^2 logN)$ まで下げることができます。

## 考え方
フィボナッチ数列の$N$項を$fib(N)$と表しこれを計算したいとします。
繰り返し二乗法のような考え方で第$N$項を計算します。そのためには以下の操作が必要になります

* $fib(A)$がわかっているときに$fib(A+1)$が計算できる
* $fib(A)$がわかっているときに$fib(2A)$が計算できる

## 最初の二項の和で表現する

実際に次の項や2倍した項を計算するために先に数列の表現を変えておきます。
$fib(A)$を計算する際に値ではなく $fib(A)=X_A \times fib(0) + Y_A \times fib(1)$ という形で表現します。
$fib(0)$と$fib(1)$の値は定義より明らかなので$X_A,Y_A$がわかれば$fib(A)$は容易に計算できます。
例えば

```math
\begin{align}
fib(2) &=1 \times fib(0)+1 \times fib(1) \\
fib(3) &=1 \times fib(2)+1 \times fib(1) \\ 
 &= (1 \times fib(0)+1 \times fib(1))+1 \times fib(1) \\ 
 &=1 \times fib(0)+2 \times fib(1)
\end{align}
```

です。
$fib(0)=0$なので$X$の数値は無駄に思えるかもしれませんが後で必要になってきます。
この変換で大事なのは以下の等式が満たされることです。
自分はこれが理解できず理解するのに時間がかかりました。

```math
fib(A+B)=X_A \times fib(B)+Y_A \times fib(B+1)
```
上の式で$B=A$にすれば$fib(2A)$が求まり、$B=1$とすれば$fib(A+1)$が求まります

## fib(A+1)を計算

具体的な数値では既に上で計算していますが式としてちゃんと計算します
既に$fib(A)=X_A \times fib(0) + Y_A \times fib(1)$となる$X_A,Y_A$がわかっているため以下のように計算できます。

```math
\begin{align}
fib(A+1)&=X_A \times fib(1) + Y_A \times fib(2) \\ 
&= X_A \times fib(1) + Y_A( 1 \times fib(0) + 1 \times fib(1) ) \\ 
&= (Y_A \times 1)fib(0) + (Y_A \times 1 + X_A)fib(1) 
\end{align}
```
式の上での$1$はフィボナッチ数列の場合は$1$ですが数列によっては違うので省略せずに書いています。

## fib(2A)を計算
上記の方法で$fib(A+1)$を先に計算しておきます。
そうすることで以下の式で計算できます。

```math
\begin{align}
fib(A+A)&=X_A \times fib(A) + Y_A \times fib(A+1) \\ 
&= X_A \times (X_A \times fib(0) + Y_A \times fib(1)) + Y_A((Y_A \times 1)fib(0) + (Y_A \times 1 + X_A)fib(1) ) \\ 
&= (X_A^2 +Y_A^2 \times 1)fib(0) + (2 \times X_A \times Y_A + Y_A^2 \times 1 )fib(1)
\end{align}
```
具体例として書きましたが。一行目の式から後はプログラムで計算させてください

# 検算
というわけで上記方法で$fib(3)$から$fib(4)$,$fib(3)$から$fib(6)$を計算します。
まず前提として$X_3=1,Y_3=2$です。

```math
\begin{align}
fib(3+1)&=(Y_3 \times 1)fib(0) + (Y_3 \times 1 + X_3)fib(1)  \\ 
&= 2 \times fib(0) + 3 \times fib(1) \\ 
&= 3
\end{align}
```

```math
\begin{align}
fib(3+3)&=(X_3^2 +Y_3^2 \times 1)fib(0) + (2 \times X_3 \times Y_3 + Y_3^2 \times 1 )fib(1) \\ 
&= 5 \times fib(0) + 8 \times fib(1) \\ 
&= 8
\end{align}
```
正しそうです
後は計算したい$N$について大きいbitから見ていってbitが立っているときに$+1$と左シフトで$\times 2$を用いてNになるように計算していけば計算ができます。


# 一般化
今までフィボナッチ数列という前提で計算してきましたが実際にはもう少し広い数列上で計算できます。
具体的には以下の線形漸化式で表現できるものが計算できます。

```math
a_{i+k} = d_0 \times a_i + d_1 \times a_{i+1}+... +d_{k-1} \times a_{i+k-1}
```
これは、フィボナッチ数列では$a_{A-1} \times 1 + a_{A-2} \times 1$と前の項と前の前の項を1倍したものを足していましたがk項前までの足し算になる、1倍とは限らないというような汎用的な表現に変わっています
伴いフィボナッチ数列では$fib(A)=X_A \times fib(0) + Y_A \times fib(1)$と表現していたものを

```math
f(A)=\{X_0,X_1,...X_{k-1}\}となる関数fが存在し、その時 \\ 
a_A = X_0 \times a_0 + .... + X_{k-1} \times a_{k-1} を満たす
```
というような感じで表現します$X_A,Y_A$だったのがk個の数値に変わった感じですね
考え方は同じで計算がややこしくなっただけです。
$f(A+1)$は以下のように計算します。

```math
\begin{align}
a_{n+1} &= X_0 \times a_1 + X_1 \times a_2 + ... + X_{k-1} \times a_k \\
&= X_0 \times a_1 + X_1 \times a_2 + ... + X_{k-1} \times (d_0 \times a_0 + d_1 \times a_1 + ... + d_{k-1} \times a_{k-1}) \\
&= (X_{k-1} \times d_0)a_0 + (X_0 + X_{k-1} \times d_1)a_1 + ... + (X_{k-1} \times d_{k-1})a_{k-1}
\end{align}
```

$f(2A)$ は以下です$a_{n+1},a_{n+2},...,a_{n+k-1}$は$f(A+1)$を求める方法で求めておきます

```math
\begin{align}
a_{n+n} &= X_0 \times a_n + X_1 \times a_n+1 + ... + X_{k-1} \times a_{n+k-1} \\
\end{align}
```
上記の式の$a_{n+1},a_{n+2},...,a_{n+k-1}$を$a_{1},a_{2},...,a_{k-1}$の式になるように展開して整理すればよいです。

これで一般の形式においても計算することができました。

# 実装
pythonコードです。$10^9+7$でmodをとっています。

```python
class kitamasa:
    def __init__(self, c0, d, mod):
        self.c0 = c0
        self.k = len(c0)
        self.d = d
        self.mod = mod

    def next(self, c):
        next_c = [0] * self.k
        next_c[0] = self.d[0] * c[-1]
        for i in range(1, self.k):
            next_c[i] = c[i - 1] + self.d[i] * c[-1]
            next_c[i] %= self.mod
        next_c[0] %= self.mod
        return next_c

    def double(self, c):
        double_c = [0] * self.k
        cs = [c]
        for i in range(self.k - 1):
            cs += (self.next(cs[-1]),)
        for i, ci in enumerate(cs):
            for j in range(self.k):
                double_c[j] += ci[j] * c[i]
                double_c[j] %= self.mod
        return double_c

    def calc(self, n):
        c = [0]*self.k
        c[1]=1
        msb = 0
        for i in range(64, -1, -1):
            if n & (1 << i):
                msb = i
                break
        for i in range(msb - 1, -1, -1):
            c = self.double(c)
            if n & (1 << i):
                c = self.next(c)
        return sum(c[i]*self.d[i] for i in range(self.k))%self.mod



# fib
k = kitamasa([0, 1], [1, 1], 10 ** 9 + 7)
# fib(21)を求める
print(k.calc(21-1))
```

# 終わりに
ここまで読んでくださりありがとうございます。
間違いなどありましたらご指摘よろしくお願いします。