【Java】フィボナッチ数(数列)を求める

　「エラトステネスの篩(ふるい)」同様、数列的なテーマとしてもよく取り上げられる「フィボナッチ数」。「黄金比(黄金数)」とも深い関係があり、実は日常的にもよく使われているとも言われている（自然現象からデザイン的なものまで）。

　まぁ、それはさておき(笑)、プログラミングで計算するのは意外と簡単だったりする。せっかくなので、いくつかの方法で求めてみよう。

■動的計画法
■再帰（メモ化再帰）
■ビネの公式
■任意精度を使う

■動的計画法でフィボナッチ数列を求める

　まずは漸化式を確認してみよう。Wikipedia にも載っているが、もう少しプログラミング用に変形した方がわかりやすいかも知れない。

(漸化式)
F₀ = 0,
F₁ = 1,
F_n = F_n-1 + F_n-2 (n ≧ 0)

●動的計画法でフィボナッチ数列を作る

import java.util.Arrays;

int n = 20;    //n 上限

//動的計画法
long[] dp = new long[n + 1];

dp[0] = 0;
dp[1] = 1;
for (int i = 2; i <= n; i++) {
    dp[i] = dp[i - 1] + dp[i - 2];
}

//F(n)を表示
System.out.println("F(" + n + ") = " + dp[n]);

//数列を表示
System.out.println(Arrays.toString(dp));

F(20) = 6765
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765]

　一般的には自然数（0 は含まない整数）で書かれているものも多いが、コンピュータで数列を作る場合は 0 も入れておいた方が何かと便利だろう。インデクスに対応できる利点もある。漸化式をそのままコーディングした感じになるので、説明はいらないだろう。

　気をつけないといけないのは、フィボナッチ数は非常に大きな値となるので、long 型で実装した場合、n = 93 でオーバーフローすることだ。それ以上の値の場合は、後述する BigDecimal（任意精度実数）を使うのも良いだろう。ちなみに、BigInteger（任意精度整数）でも構わないが、BigDecimal の方がわずかに速いようだ。


■再帰でフィボナッチ数(数列)を求める

　次に再帰で フィボナッチ数 F(n) を求めてみよう。再帰でもとてもシンプルに書けることがわかる。

●再帰でフィボナッチ数を求める

int n = 20;

//F(n)を表示
System.out.println("F(" + n + ") = " + fibonacci(n));

//再帰関数
long fibonacci(int n) {
    if (n <= 1) {
        return  n;
    }
    return fibonacci(n - 1) + fibonacci(n - 2);
}

F(20) = 6765

　ただ、上記の方法だと n = 40 ほどになると非常に重くなる。なぜなら、再帰で足しあわせている部分「fibonacci(n - 1) + fibonacci(n - 2)」は入れ子になっていて、同じ計算を何度も繰り返しているため、無駄が多いからだ。この関数は引数が同じなら戻値も同じになるので、メモ化することで高速化を図れる。またメモは数列にもなるので表示してみよう。

●メモ化再帰でフィボナッチ数列を作る

import java.util.Arrays;

int n = 20;

long[] fib = new long[n + 1];  //メモ(グローバルに置く)

//F(n)を表示
System.out.println("F(" + n + ") = " + fibonacci(n));

//数列を表示
System.out.println(Arrays.toString(fib));

//メモ化再帰
long fibonacci(int n) {
    if (n <= 1) {
        return fib[n] = n;
    }
    if (fib[n] > 0) {
        return fib[n];
    }
    return fib[n] = fibonacci(n - 1) + fibonacci(n - 2);
}

F(20) = 6765
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765]

　この方法なら n = 40 以上でも一瞬で求められる。

　ただし、上記の動的計画法でも再帰でも、long で実装した場合、 n = 93 でオーバーフローになるので気をつけよう。


■ビネの公式でフィボナッチ数を求める

　もう１つ、一般項を求める式（ビネの公式）を用いてフィボナッチ数を求めてみよう。

●ビネの公式でフィボナッチ数を求める

final double GOLDEN_RATIO = 1.618033988749895;  //黄金比
final double SQRT5 = Math.sqrt(5);  //√5

int n = 20;

//ビネの公式
//(上の式) n = 72 から誤差が出る
double ans1 = (Math.pow(1.0 + SQRT5, n) - Math.pow(1.0 - SQRT5, n)) / (Math.pow(2, n) * SQRT5);
System.out.println("F(" + n + ") = " + (long)ans1 + " (" + ans1 + ")");

//(下の式) n = 71 から誤差が出る
double ans2 = Math.floor(Math.pow(GOLDEN_RATIO, n) / SQRT5 + 0.5);
System.out.println("F(" + n + ") = " + (long)ans2 + " (" + ans2 + ")");

F(20) = 6765 (6765.000000000005)
F(20) = 6765 (6765.0)

　ただし、これらは近似式のため誤差が出るので注意しよう。 double で実装した場合、上の式では n = 72 から、下の式では n = 71 から誤差が出る。式の解説は Wikipedia に載っているので参照して欲しいが、プログラミングするなら下の方の式が使いやすい。カッコの上の部分が無いような記号は床関数の記号で、小数点切り捨てと考えれば、Math.floor() でも整数型のキャストでも良いだろう（正の整数の場合）。よくわからない記号は「数学記号の表」などで調べると良い。

　ビネの公式での誤差の出る範囲を抑える１つの方法は、BigDecimal のような任意精度実数を使って、桁の精度を上げることである。試しに黄金比や √5 の値を 50 桁まで上げて、先ほどの式（下の方の式、上の方の式は上手くいかない）を BigDecimal に書き換えてみると、n = 231 までは誤差が出ない。

●ビネの公式を BigDecimal で書いてみる(小数点以下 50桁)

import java.math.BigDecimal;
import java.math.RoundingMode;

int n = 231;  //232 から誤差が出る

/** 任意精度実数の黄金比 (50桁) */
final BigDecimal BIG_DECIMAL_GOLDEN = new BigDecimal("1.61803398874989484820458683436563811772030917980576");

/** 任意精度実数の √5 (50桁) */
final BigDecimal BIG_DECIMAL_SQRT5 = new BigDecimal("2.23606797749978969640917366873127623544061835961152");

BigDecimal ans = BIG_DECIMAL_GOLDEN.pow(n).divide(BIG_DECIMAL_SQRT5, 50, RoundingMode.FLOOR).add(new BigDecimal(0.5));

System.out.println("F(" + n + ") = " + ans.setScale(0, RoundingMode.FLOOR) + " (" + ans + ")");

F(231) = 844617150046923109759866426342507997914076076194 (844617150046923109759866426342507997914076076194.15703878686668096410036906990688047595887460281556)

　ビネの公式も黄金比や √5 の精度を上げれば、もっと大きな n でも誤差が出ないようにできると思うが、実際には計算が非常に重くなるのであまり実用的ではないかも知れない。その場合、はじめに書いた動的計画法を BigDecimal に書き換えた方が実行速度も速く、誤差もない。

●はじめの動的計画法を BigDecimal に書き換えたもの

import java.math.BigDecimal;
import java.util.Arrays;

int n = 100;

BigDecimal[] dp = new BigDecimal[n + 1];

dp[0] = BigDecimal.ZERO;
dp[1] = BigDecimal.ONE;
for (int i = 2; i <= n; i++) {
    dp[i] = dp[i - 1].add(dp[i - 2]);
}

//F(n)を表示
System.out.println("F(" + n + ") = " + dp[n]);

//数列を表示
System.out.println(Arrays.toString(dp));

F(100) = 354224848179261915075
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418, 317811, 514229, 832040, 1346269, 2178309, 3524578, 5702887, 9227465, 14930352, 24157817, 39088169, 63245986, 102334155, 165580141, 267914296, 433494437, 701408733, 1134903170, 1836311903, 2971215073, 4807526976, 7778742049, 12586269025, 20365011074, 32951280099, 53316291173, 86267571272, 139583862445, 225851433717, 365435296162, 591286729879, 956722026041, 1548008755920, 2504730781961, 4052739537881, 6557470319842, 10610209857723, 17167680177565, 27777890035288, 44945570212853, 72723460248141, 117669030460994, 190392490709135, 308061521170129, 498454011879264, 806515533049393, 1304969544928657, 2111485077978050, 3416454622906707, 5527939700884757, 8944394323791464, 14472334024676221, 23416728348467685, 37889062373143906, 61305790721611591, 99194853094755497, 160500643816367088, 259695496911122585, 420196140727489673, 679891637638612258, 1100087778366101931, 1779979416004714189, 2880067194370816120, 4660046610375530309, 7540113804746346429, 12200160415121876738, 19740274219868223167, 31940434634990099905, 51680708854858323072, 83621143489848422977, 135301852344706746049, 218922995834555169026, 354224848179261915075]

　最後に long, double (ビネの公式：下の式), BigDecimal で計算した結果をダンプしてみよう。オーバーフローや誤差を確認しておけば、どれを使うかの目安になるかも知れない。値を確認するには以下のようなサイトを参考にすれば良い。

（参考）[数列]
・１００番目までのフィボナッチ数列
・１番目から１０００番目までのフィボナッチ数列

F(0)
long       = 0
double     = 0
BigDecimal = 0
F(1)
long       = 1
double     = 1
BigDecimal = 1
F(2)
long       = 1
double     = 1
BigDecimal = 1
F(3)
long       = 2
double     = 2
BigDecimal = 2
F(4)
long       = 3
double     = 3
BigDecimal = 3
F(5)
long       = 5
double     = 5
BigDecimal = 5
F(6)
long       = 8
double     = 8
BigDecimal = 8
F(7)
long       = 13
double     = 13
BigDecimal = 13
F(8)
long       = 21
double     = 21
BigDecimal = 21
F(9)
long       = 34
double     = 34
BigDecimal = 34
F(10)
long       = 55
double     = 55
BigDecimal = 55
F(11)
long       = 89
double     = 89
BigDecimal = 89
F(12)
long       = 144
double     = 144
BigDecimal = 144
F(13)
long       = 233
double     = 233
BigDecimal = 233
F(14)
long       = 377
double     = 377
BigDecimal = 377
F(15)
long       = 610
double     = 610
BigDecimal = 610
F(16)
long       = 987
double     = 987
BigDecimal = 987
F(17)
long       = 1597
double     = 1597
BigDecimal = 1597
F(18)
long       = 2584
double     = 2584
BigDecimal = 2584
F(19)
long       = 4181
double     = 4181
BigDecimal = 4181
F(20)
long       = 6765
double     = 6765
BigDecimal = 6765
F(21)
long       = 10946
double     = 10946
BigDecimal = 10946
F(22)
long       = 17711
double     = 17711
BigDecimal = 17711
F(23)
long       = 28657
double     = 28657
BigDecimal = 28657
F(24)
long       = 46368
double     = 46368
BigDecimal = 46368
F(25)
long       = 75025
double     = 75025
BigDecimal = 75025
F(26)
long       = 121393
double     = 121393
BigDecimal = 121393
F(27)
long       = 196418
double     = 196418
BigDecimal = 196418
F(28)
long       = 317811
double     = 317811
BigDecimal = 317811
F(29)
long       = 514229
double     = 514229
BigDecimal = 514229
F(30)
long       = 832040
double     = 832040
BigDecimal = 832040
F(31)
long       = 1346269
double     = 1346269
BigDecimal = 1346269
F(32)
long       = 2178309
double     = 2178309
BigDecimal = 2178309
F(33)
long       = 3524578
double     = 3524578
BigDecimal = 3524578
F(34)
long       = 5702887
double     = 5702887
BigDecimal = 5702887
F(35)
long       = 9227465
double     = 9227465
BigDecimal = 9227465
F(36)
long       = 14930352
double     = 14930352
BigDecimal = 14930352
F(37)
long       = 24157817
double     = 24157817
BigDecimal = 24157817
F(38)
long       = 39088169
double     = 39088169
BigDecimal = 39088169
F(39)
long       = 63245986
double     = 63245986
BigDecimal = 63245986
F(40)
long       = 102334155
double     = 102334155
BigDecimal = 102334155
F(41)
long       = 165580141
double     = 165580141
BigDecimal = 165580141
F(42)
long       = 267914296
double     = 267914296
BigDecimal = 267914296
F(43)
long       = 433494437
double     = 433494437
BigDecimal = 433494437
F(44)
long       = 701408733
double     = 701408733
BigDecimal = 701408733
F(45)
long       = 1134903170
double     = 1134903170
BigDecimal = 1134903170
F(46)
long       = 1836311903
double     = 1836311903
BigDecimal = 1836311903
F(47)
long       = 2971215073
double     = 2971215073
BigDecimal = 2971215073
F(48)
long       = 4807526976
double     = 4807526976
BigDecimal = 4807526976
F(49)
long       = 7778742049
double     = 7778742049
BigDecimal = 7778742049
F(50)
long       = 12586269025
double     = 12586269025
BigDecimal = 12586269025
F(51)
long       = 20365011074
double     = 20365011074
BigDecimal = 20365011074
F(52)
long       = 32951280099
double     = 32951280099
BigDecimal = 32951280099
F(53)
long       = 53316291173
double     = 53316291173
BigDecimal = 53316291173
F(54)
long       = 86267571272
double     = 86267571272
BigDecimal = 86267571272
F(55)
long       = 139583862445
double     = 139583862445
BigDecimal = 139583862445
F(56)
long       = 225851433717
double     = 225851433717
BigDecimal = 225851433717
F(57)
long       = 365435296162
double     = 365435296162
BigDecimal = 365435296162
F(58)
long       = 591286729879
double     = 591286729879
BigDecimal = 591286729879
F(59)
long       = 956722026041
double     = 956722026041
BigDecimal = 956722026041
F(60)
long       = 1548008755920
double     = 1548008755920
BigDecimal = 1548008755920
F(61)
long       = 2504730781961
double     = 2504730781961
BigDecimal = 2504730781961
F(62)
long       = 4052739537881
double     = 4052739537881
BigDecimal = 4052739537881
F(63)
long       = 6557470319842
double     = 6557470319842
BigDecimal = 6557470319842
F(64)
long       = 10610209857723
double     = 10610209857723
BigDecimal = 10610209857723
F(65)
long       = 17167680177565
double     = 17167680177565
BigDecimal = 17167680177565
F(66)
long       = 27777890035288
double     = 27777890035288
BigDecimal = 27777890035288
F(67)
long       = 44945570212853
double     = 44945570212853
BigDecimal = 44945570212853
F(68)
long       = 72723460248141
double     = 72723460248141
BigDecimal = 72723460248141
F(69)
long       = 117669030460994
double     = 117669030460994
BigDecimal = 117669030460994
F(70)
long       = 190392490709135
double     = 190392490709135
BigDecimal = 190392490709135
F(71)
long       = 308061521170129
double     = 308061521170130      //← ここから誤差が出る
BigDecimal = 308061521170129
F(72)
long       = 498454011879264
double     = 498454011879265
BigDecimal = 498454011879264
F(73)
long       = 806515533049393
double     = 806515533049395
BigDecimal = 806515533049393
F(74)
long       = 1304969544928657
double     = 1304969544928660
BigDecimal = 1304969544928657
F(75)
long       = 2111485077978050
double     = 2111485077978055
BigDecimal = 2111485077978050
F(76)
long       = 3416454622906707
double     = 3416454622906716
BigDecimal = 3416454622906707
F(77)
long       = 5527939700884757
double     = 5527939700884772
BigDecimal = 5527939700884757
F(78)
long       = 8944394323791464
double     = 8944394323791488
BigDecimal = 8944394323791464
F(79)
long       = 14472334024676221
double     = 14472334024676260
BigDecimal = 14472334024676221
F(80)
long       = 23416728348467685
double     = 23416728348467744
BigDecimal = 23416728348467685
F(81)
long       = 37889062373143906
double     = 37889062373144008
BigDecimal = 37889062373143906
F(82)
long       = 61305790721611591
double     = 61305790721611752
BigDecimal = 61305790721611591
F(83)
long       = 99194853094755497
double     = 99194853094755776
BigDecimal = 99194853094755497
F(84)
long       = 160500643816367088
double     = 160500643816367552
BigDecimal = 160500643816367088
F(85)
long       = 259695496911122585
double     = 259695496911123328
BigDecimal = 259695496911122585
F(86)
long       = 420196140727489673
double     = 420196140727490880
BigDecimal = 420196140727489673
F(87)
long       = 679891637638612258
double     = 679891637638614270
BigDecimal = 679891637638612258
F(88)
long       = 1100087778366101931
double     = 1100087778366105090
BigDecimal = 1100087778366101931
F(89)
long       = 1779979416004714189
double     = 1779979416004719360
BigDecimal = 1779979416004714189
F(90)
long       = 2880067194370816120
double     = 2880067194370824700
BigDecimal = 2880067194370816120
F(91)
long       = 4660046610375530309
double     = 4660046610375544800
BigDecimal = 4660046610375530309
F(92)
long       = 7540113804746346429
double     = 7540113804746369000
BigDecimal = 7540113804746346429
F(93)
long       = -6246583658587674878      //← ここでオーバーフロー
double     = 12200160415121914000
BigDecimal = 12200160415121876738
F(94)
long       = 1293530146158671551
double     = 19740274219868283000
BigDecimal = 19740274219868223167
F(95)
long       = -4953053512429003327
double     = 31940434634990200000
BigDecimal = 31940434634990099905
F(96)
long       = -3659523366270331776
double     = 51680708854858490000
BigDecimal = 51680708854858323072
F(97)
long       = -8612576878699335103
double     = 83621143489848690000
BigDecimal = 83621143489848422977
F(98)
long       = 6174643828739884737
double     = 135301852344707190000
BigDecimal = 135301852344706746049
F(99)
long       = -2437933049959450366
double     = 218922995834555900000
BigDecimal = 218922995834555169026
F(100)
long       = 3736710778780434371
double     = 354224848179263100000
BigDecimal = 354224848179261915075

（関連記事）
【Java】任意精度整数(多倍長整数)で演算する

■参考になる書籍