今日は12月3日。
1・2・3と数字が並んでおり、フィボナッチ数列みたいといえる。
そこで、今回は IchigoJam でフィボナッチ数列の計算をしてみる。

※IchigoJamはjig.jpの登録商標です。

フィボナッチ数列

(この記事における) フィボナッチ数列は

最初の2項は $a_0 = 0, a_1 = 1$
その後の項は $a_n = a_{n-1} + a_{n-2}$ ($n \geq 2$)

で表される数列である。
すなわち、$a_2$ 以降はその前の2項の和である。

素直なプログラムで求める

たとえば、以下のプログラムでフィボナッチ数列の N 項目 $a_N$ を十進数で表したときの下4桁を求めることができる。
(値はどんどん大きくなるため、今回は下4桁のみを求める)

ループを用いて、前の項 P (prev) と今の項 C (current) の値を更新していく。

10 ' ﾌｨﾎﾞﾅｯﾁ ｽｳﾚﾂ (ｸﾞﾁｮｸ)
20 INPUT N
30 CLT
40 P=1:C=0:IF N=0 GOTO 80
50 FOR I=1 TO N
60 V=P+C:P=C:C=V%10000
70 NEXT
80 ?"N=";N;",RES=";C;",TIME=";TICK()

行列を用いて求める

フィボナッチ数列を行列で表す

フィボナッチ数列の計算は、$n \geq 1$ のとき

\begin{eqnarray}
a_{n+2} &=& a_{n+1} + a_n \\
a_{n+1} &=& a_{n+1}
\end{eqnarray}

すなわち

\begin{pmatrix}
a_{n+2} \\
a_{n+1}
\end{pmatrix}
=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix}
a_{n+1} \\
a_n
\end{pmatrix}

と表せる。同様に

\begin{pmatrix}
a_{n+3} \\
a_{n+2}
\end{pmatrix}
=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix}
a_{n+2} \\
a_{n+1}
\end{pmatrix}

なので、

\begin{pmatrix}
a_{n+3} \\
a_{n+2}
\end{pmatrix}
=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}^2
\begin{pmatrix}
a_{n+1} \\
a_n
\end{pmatrix}

となる。これを繰り返すと、$a_1 = 1, a_0 = 0$ を用いて

\begin{pmatrix}
a_{n+1} \\
a_n
\end{pmatrix}
=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}^n
\begin{pmatrix}
1 \\
0
\end{pmatrix}

という関係が導ける。

累乗を効率よく計算する

前の節で求めた関係より、

\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}^n

を効率よく求めることができれば、フィボナッチ数列を効率よく計算できることがわかる。
そこで、繰り返し二乗法を用いてこれを効率よく求める。

繰り返し二乗法とは、何乗するかの値を2の累乗の和で表し、「もとの値の(2の累乗)乗」と「それらのうち使うものの積」を計算していくことで、もとの値の累乗を計算するアルゴリズムである。
具体的には、以下のように計算を行う。

擬似コード

// A：何の累乗を求めるか
// N：何乗を求めるか (非負整数)
function pow(A, N) {
  R ← 1 // 計算結果 (「1」は単位元 (何にかけてももとの値になる値) を表す)
  D ← A // もとの値の(2の累乗)乗 (次に掛ける可能性がある値)
  while (N > 0) {     // N が正の間、以下を繰り返す
    if (N % 2 == 1) { // N を 2 で割った余りが 1 ならば
      R ← R * D       // 計算結果にもとの値の(2の累乗)乗を掛ける
    }
    D ← D * D         // もとの値の(2の累乗)乗を更新する
    N ← N / 2         // N を 2 で割る (小数点以下は切り捨て)
  }
  return R // 計算結果を返す
}

N を 2 で割っていくことで 2 進数のそれぞれの桁を最下位に持っていき、それが 0 (掛けない) か 1 (掛ける) かの判定を行う。
それとともに、D を二乗していくことで、もとの値の2乗、4乗 (2乗の2乗)、8乗 (4乗の2乗)、…を計算し、掛ける準備をする。

たとえば、A = 3, N = 13 として $3^{13}$ を求める流れは、以下のようになる。

N	N (2進数)	D	R
13	1101	3	1×3=3
6	110	3×3=9	3
3	11	9×9=81	3×81=243
1	1	81×81=6561	243×6561=1594323

実装

累乗の結果を

\begin{pmatrix}
A & B \\
C & D
\end{pmatrix}

に、次に掛けるかもしれない行列を

\begin{pmatrix}
E & F \\
G & H
\end{pmatrix}

に格納する。
このとき、累乗の結果の更新は

\begin{pmatrix}
O & P \\
Q & R
\end{pmatrix}
=
\begin{pmatrix}
A & B \\
C & D
\end{pmatrix}
\times
\begin{pmatrix}
E & F \\
G & H
\end{pmatrix}
=
\begin{pmatrix}
AE+BG & AF+BH \\
CE+DG & CF+DH
\end{pmatrix}

次に掛けるかもしれない行列の更新は

\begin{pmatrix}
O & P \\
Q & R
\end{pmatrix}
=
\begin{pmatrix}
E & F \\
G & H
\end{pmatrix}
\times
\begin{pmatrix}
E & F \\
G & H
\end{pmatrix}
=
\begin{pmatrix}
E^2+FG & EF+FH \\
EG+GH & FG+H^2
\end{pmatrix}

で計算できる。
累乗の結果の初期値は単位行列

\begin{pmatrix}
A & B \\
C & D
\end{pmatrix}
=
\begin{pmatrix}
1 & 0 \\
0 & 1
\end{pmatrix}

であり、次に掛けるかもしれない行列の初期値は関係式に出てきた行列

\begin{pmatrix}
E & F \\
G & H
\end{pmatrix}
=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}

である。

また、掛け算をする際、十進数で4桁の整数同士をそのまま掛けると、32,767 までの整数しか扱えない IchigoJam BASIC ではオーバーフローする可能性がある。
そこで、繰り返し二乗法と同じ要領で、掛け算のかわりに足し算、二乗のかわりに二倍をすることで、1ビットずつ筆算のように掛け算を行う。
すると、最大で 10,000 の2倍の 20,000 程度の数までしか扱わずにすみ、IchigoJam BASIC で計算ができるようになる。

10 ' ﾌｨﾎﾞﾅｯﾁ ｽｳﾚﾂ (ｷﾞｮｳﾚﾂ)
20 INPUT N
30 CLT
40 M=10000:GOTO 80
50 Z=0
60 IF Y&1 Z=(Z+X)%M
70 X=X*2%M:Y=Y/2:IF Y GOTO 60 ELSE RETURN
80 A=1:B=0:C=0:D=1
90 E=1:F=1:G=1:H=0
100 K=N
110 IF K&1=0 GOTO 170
120 X=A:Y=E:GOSUB50:O=Z:X=B:Y=G:GOSUB50:O=O+Z
130 X=A:Y=F:GOSUB50:P=Z:X=B:Y=H:GOSUB50:P=P+Z
140 X=C:Y=E:GOSUB50:Q=Z:X=D:Y=G:GOSUB50:Q=Q+Z
150 X=C:Y=F:GOSUB50:R=Z:X=D:Y=H:GOSUB50:R=R+Z
160 A=O%M:B=P%M:C=Q%M:D=R%M
170 X=E:Y=E:GOSUB50:O=Z:X=F:Y=G:GOSUB50:O=O+Z
180 X=E:Y=F:GOSUB50:P=Z:X=F:Y=H:GOSUB50:P=P+Z
190 X=E:Y=G:GOSUB50:Q=Z:X=G:Y=H:GOSUB50:Q=Q+Z
200 X=F:Y=G:GOSUB50:R=Z:X=H:Y=H:GOSUB50:R=R+Z
210 E=O%M:F=P%M:G=Q%M:H=R%M
220 K=K/2:IF K GOTO 110
230 ?"N=";N;",RES=";C;",TIME=";TICK()

50～70行目で、乗算結果の下4桁を求めるサブルーチンを定義している。
このサブルーチンは、X×Y の下4桁を求め、Z に格納する。
計算の際、X および Y の値は破壊される。

80～100行目で、変数の初期値を設定する。
110行目で、繰り返し二乗法の乗算を行うかを判定する。
120～160行目で、繰り返し二乗法の乗算 (累乗の結果の更新) を行う。
170～210行目で、繰り返し二乗法の二乗 (次に掛けるかもしれない行列の更新) を行う。
220行目で、繰り返し二乗法の状態を更新し、ループを進める。

計算時間の比較

IchigoJam Q (IchigoJam BASIC 1.4.3) でそれぞれのプログラムを実行し、計算時間を比較した。
そのために、それぞれのプログラムに

20 FOR N=1000 TO 32000 STEP 1000
300 NEXT

を追加し、N の値を自動で設定して計算できるようにした。
また、事前に VIDEO 0 を実行し、映像出力を切った状態で1回測定を行った。
結果は以下のようになった。

詳細

N の値、およびそのときの各プログラムによる計算時間 (TICK() の単位) を以下に示す。

N	素直なプログラム	行列を用いたプログラム
1000	149	276
2000	299	307
3000	449	337
4000	597	332
5000	747	319
6000	896	375
7000	1046	359
8000	1195	354
9000	1344	342
10000	1494	350
11000	1643	408
12000	1792	403
13000	1942	367
14000	2091	399
15000	2240	384
16000	2390	380
17000	2539	369
18000	2688	377
19000	2838	387
20000	2986	381
21000	3136	345
22000	3285	444
23000	3435	429
24000	3584	424
25000	3734	388
26000	3883	397
27000	4032	429
28000	4181	424
29000	4330	388
30000	4480	421
31000	4630	404
32000	4779	399

素直に計算した場合、N の値に比例する感じで、どんどん計算時間が伸びた。
一方、行列を用いて計算した場合、N の値が小さいときは素直に計算するよりも計算時間が長くなったものの、N の値が大きくなっても計算時間があまり伸びず、N の値によっては素直に計算した場合の10分の1以下の時間で計算を行うことができた。

ちなみに

計算してみると

\begin{eqnarray}
a_{15000} \bmod 10000 &=& 0 \\
a_{15001} \bmod 10000 &=& 1
\end{eqnarray}

となるので、出現する値がループし、「フィボナッチ数列の十進数で下4桁」は、第何項かの値を 15000 で割ってから求めても同じ値になる。
この周期は、他の値で余りを取った場合は変わることがあるだろう。