動的計画法 - ヽ｜∵｜ゝ(Fantom) の開発blog？

このページの記事一覧

2016/07/10 【Java】フィボナッチ数(数列)を求める
2016/02/01 【Java】パスカルの三角形を使って組合せを求める
2015/04/28 【Java】ナップサック問題(knapsack)[動的計画法]

« prev next »

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

2016/07/10 Sun [edit]

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

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

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

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

　まずは漸化式を確認してみよう。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】任意精度整数(多倍長整数)で演算する

■参考になる書籍

関連記事

スポンサーサイト

category: Java

thread: プログラミング

janre: コンピュータ

tag: 算術関数数列動的計画法検証

tb: 0 cm: --

【Java】パスカルの三角形を使って組合せを求める

2016/02/01 Mon [edit]

　前回、漸化式を使っての組合せの数を求めたが、別の解法として「パスカルの三角形」を利用して組合せを求める方法もある。パスカルの三角形とは「二項係数」とか「二項定理」とも呼ばれ、高校数学でやった二項式を展開するときの係数に合致しているが、他にも色々利用できる性質がある。

　まずは簡単に作り方を覚えてしまおう。非常に単純なので一度聞いたら忘れることもないだろう。作り方は以下のようになる。

１．頂点（一番上）に「１」を書き、ここをスタート地点とする。

２．スタート地点から左右（矢印の方向：斜め下）に向かって進んでいく。その際、交差する地点はその１つ前の値を合計する。これを繰り返す。

　これは道に例えると「網の目のように必ず２つの分岐点がある道を進むとき、その地点に最短で辿り着く方法は何通りあるか？」という数になる。つまり、地点Aに辿り着く方法が２通り、地点Bに辿り着く方法が１通りあるとき、その地点AとBが合流する地点Cは 2+1=3 通りあるという単純なものである。

　そしてこの値は、以下のように組合せを求める式の答えとも合致している。

　つまりこのパスカルの三角形を作ることをシミュレーションすれば、組合せの数も簡単に求めることができるというわけだ。さっそく上記の作り方をそのままコーディングしてみよう。

●パスカルの三角形を作る [O(hw)]

/**<h1>パスカルの三角形を作る</h1>
 * <p>パスカルの三角形を斜めに傾けた配列を作成。左上が(0,0)で右下に拡がる感じになる。</p>
 * @param w ： 横数 (1～33)
 * @param h ： 縦数 (1～33)
 * @return<b>long[][]</b> ： パスカルの三角形を配列化したもの
 */
public static final long[][] pascalsTriangle(final int w, final int h) {
    final long[][] arr = new long[h + 1][w + 1];
    arr[0][0] = 1;

    for (int i = 0; i <= h; i++) {
        for (int j = 0; j <= w; j++) {
            if (j > 0) {
                arr[i][j] += arr[i][j - 1];
            }
            if (i > 0) {
                arr[i][j] += arr[i - 1][j];
            }
        }
    }

    return arr;
}


//メインでは...
int w = 5;
int h = 4;
long[][] pascal = pascalsTriangle(w, h);    //パスカルの三角形

　パスカルの三角形を下の図のように斜め45度傾けて配列に入れたと考えれば簡単だろう。

　つまり組合せの数を求めたいなら、この表を作って、値を取得すれば良いということになる。例えば、5C3 を求めたければ、配列の pascal[3][2] = 10 が「5個のものから3個取り出す組合せ」の数になる。nCr に置き換えれば、pascal[n - r][r] となる。また表を見てわかるように、対角線で値が対称になっているので、5C3 = 5C2 (nCr = nCn-r) となり、配列の縦横を入れ替えても同じとなる。もし、関数として作って置きたいなら、以下のような感じにすれば良い。

●パスカルの三角形から組合せの数を求める [O(n²)]

/**<h1>n 個のものから r 個取り出した組合せの数を求める</h1>
 * <p>nCr をパスカルの三角形から求める。
 * <br>n, r が 大きすぎるとオーバーフローするので注意。</p>
 * @param n ： すべての数 (1～66) (負の値はエラー)
 * @param r ： 取り出す数 (1～33) (負の値はエラー)
 * @return<b>long</b> ： 組合せの数
 */
public static final long combination(final int n, final int r) {
    if (n < 0 || r < 0 || n < r) {
        throw new ArithmeticException("n = " + n + ", r = " + r);
    }
    if (r == 0 || r == n) {
        return 1;
    }
    if (r == 1) {
        return n;
    }

    //パスカルの三角形
    final long[][] pascal = pascalsTriangle(r, n - r);
    return pascal[n - r][r];
}

　漸化式での組合せ[O(n)]と比べたメリットは、計算が単純な加算しかしていないため、オーバフローが少しばかりしにくい所だろうか。また、大量に計算するなら、１つパスカルの三角形を作っておいて[O(n²)]、配列から値を取り出すだけにすれば[O(1)]、実行負荷はほとんど無くなる（例えば、30個の組合せを1000回計算するなら、漸化式では 30[O(n)]x1000回 = 30000回演算、パスカルの三角形なら 30x30[O(n²)]=900 +1000[O(1)]回 = 1900回演算で、計算量が少なくなるので速い）。しかし単体で計算するなら、毎回配列を作る[O(n²)]ので、実行速度は遅くなるので注意しよう。

　ここで試しに練習問題をやってみよう。これはどちらかというとプログラミング問題というより数学の問題だ。

４人でじゃんけんを１回するとき、あいこになる確率はどれくらいか？

//パスカルの三角形での解法
int n = 4;    //4人
long[][] pascal = pascalsTriangle(n, n); //大きさは適当(n x n あればカバーできる)
long sum = 0;    //勝敗がつく事象の合計
for (int i = 1; i < n; i++) {    //4C1～4C3 まで(4人のうち1～3人勝つ[=勝敗がつく]組合せ)
    sum += pascal[n - i][i];
}
sum *= 3;    //グー、チョキ、パーの３手でそれぞれ勝てるため、３倍する
double d = 1.0 - sum / Math.pow(3, n);    //あいこ = 全事象 - 勝敗がつく事象 [=余事象]
System.out.println(d);    //0.4814814814814815 (13/27)

　答えは約48% (= 13/27) である。「勝ち・負け・あいこの３つだから、1/3 じゃないの？」と思った人は気をつけよう。それは前回の「コイン４枚のうち２枚が表の確率」同様、ただ単に思い込みをしているだけだ。

　簡単に解法を書いてみると、じゃんけんの手はグー・チョキ・パーの３手だから、すべての事象は 3 x 3 x 3 x 3 = 81 通りである。このうち勝敗がつくのは４人のうち１～３人の誰かが勝てばいいわけだから、4C1 + 4C2 + 4C3 の合計になり、それぞれ３手あるので３倍になる。「あいこ」になる確率は「勝敗がつかない」確率だから、余事象を用いて (1 - 勝敗がつく事象) で、1 - 14/27 = 13/27 ≒ 0.48 (48%) となる。

　ちなみに、この問題はわざわざパスカルの三角形を使った解法で試してみたが、「n人でじゃんけんを１回するとき、あいこになる確率」を一般化すると、

となるので、一発計算で求めることもできる。

//あいこになる確率を一般化したものでの解法
int n = 4;    //4人
double d = 1.0 - (Math.pow(2, n) - 2.0) / Math.pow(3, n - 1);    //あいこ = 全事象 - 勝敗がつく事象 [=余事象]
System.out.println(d);    //0.4814814814814815 (13/27)

　試しに 2人～20人でじゃんけんをしたとき、あいこになる確率を計算してみよう。

2人でじゃんけんを１回するとき、あいこになる確率 = 0.33333333333333337
3人でじゃんけんを１回するとき、あいこになる確率 = 0.33333333333333337
4人でじゃんけんを１回するとき、あいこになる確率 = 0.4814814814814815
5人でじゃんけんを１回するとき、あいこになる確率 = 0.6296296296296297
6人でじゃんけんを１回するとき、あいこになる確率 = 0.7448559670781894
7人でじゃんけんを１回するとき、あいこになる確率 = 0.8271604938271605
8人でじゃんけんを１回するとき、あいこになる確率 = 0.8838591678097851
9人でじゃんけんを１回するとき、あいこになる確率 = 0.922267946959305
10人でじゃんけんを１回するとき、あいこになる確率 = 0.9480770207793527
11人でじゃんけんを１回するとき、あいこになる確率 = 0.9653508103439516
12人でじゃんけんを１回するとき、あいこになる確率 = 0.9768892501707621
13人でじゃんけんを１回するとき、あいこになる確率 = 0.9845890700943284
14人でじゃんけんを１回するとき、あいこになる確率 = 0.9897247922786035
15人でじゃんけんを１回するとき、あいこになる確率 = 0.9931494433687528
16人でじゃんけんを１回するとき、あいこになる確率 = 0.9954328228623964
17人でじゃんけんを１回するとき、あいこになる確率 = 0.9969551687804513
18人でじゃんけんを１回するとき、あいこになる確率 = 0.9979700970332521
19人でじゃんけんを１回するとき、あいこになる確率 = 0.9986467261931519
20人でじゃんけんを１回するとき、あいこになる確率 = 0.9990978157413181

　つまり人数が増えていくたびに、あいこになる確率は限りなく 100% に近づいていく。5人以上のとき 60% を超え、9人以上のときは 90% を超えるので、もはや勝負がつかなくなる。人数が増えたとき、あいこが延々と続くようになるのはこういうことだ(笑)。人数が多いとき、4人以下で分割してじゃんけんした方が効率が良いことも理にかなっている。

　また、パスカルの三角形の作り方のコードを見て、アルゴリズマーな人なら(笑)、「あれ？このコード見たことあるな…」と思ったかもしれない。実はこれは「格子状の道の経路数を求める」(ヒョウ本 P.176～177)に載っている「動的計画法」のコードと同じだったりする。　

　その問題を挙げてみると、

以下のような格子状の道があるとき、S から G まで最短で移動する方法は何通りあるか？

というような問いである。

　実はこれも高校数学の教科書に載っているのだが、解法を簡単に書いてみると、下へ進む道を a 、右に進む道を b として１列に並べると、aaaabbbbb のようになり「9本の道から4本の道を選ぶ組合せの数」になることがわかる（9本の道から4本の道を a とすれば、残りの5本の道は b と決まるため）。そして、この格子状の経路は先のパスカルの三角形を斜めに倒した表と一致するので、全く同じコードになるわけだ。

　つまり、この格子状の経路を数え上げる問題は組合せの問題と本質的に同じになる。だから漸化式での方法でも一発で求めることができることになる。

　しかし、例えばこの格子状の道に通行止めがあった場合を考えよう。そうなると単純な組合せの計算一発では求められなくなる。その場合は、パスカルの三角形の作り方を改造すれば意外と簡単に求められるのだ。例えば先のパスカルの三角形を作るコードに通行止めのフラグでも入れて「その地点には行けない(経路数 = 0)」を挿入してしまえば良い。

●パスカルの三角形を改造して経路数を求める

//パスカルの三角形に通行止めフラグ(closed[][])を入れた
public static final long[][] pascalsTriangle(final int w, final int h, final boolean[][] closed) {
    final long[][] arr = new long[h + 1][w + 1];
    arr[0][0] = 1;

    for (int i = 0; i <= h; i++) {
        for (int j = 0; j <= w; j++) {
            if (closed[i][j]) {
                continue;
            }
            if (j > 0) {
                arr[i][j] += arr[i][j - 1];
            }
            if (i > 0) {
                arr[i][j] += arr[i - 1][j];
            }
        }
    }

    return arr;
}


//メインでは・・・
int w = 5;
int h = 4;
boolean[][] closed = new boolean[h + 1][w + 1];
closed[1][3] = true;    //通行止め
closed[2][1] = true;    //通行止め
closed[3][2] = true;    //通行止め
long[][] pascal = pascalsTriangle(w, h, closed);
System.out.println(pascal[h][w]);    //25 となる

　この手の経路数数え上げの問題は再帰を使うのも常套だが、平面を再帰で探索すると、O(2^hw) くらいの計算量はかかってしまうので、非常に重くなる。あと、マスが 100 x 100 = 10000 を超えると Stack Overflow を起こすこともあるので気をつけよう。このパスカルの三角形を改造した解法（動的計画法）なら、計算量は O(hw) でしかない。

　組合せの数と経路数の数え上げのロジックが繋がっていることを覚えておくと、利用価値があるかもしれない。色々考察してみると良いだろう。

（関連記事）
【Java】組合せの数を求める
【Java】順列の数を求める
【Java】プリミティブ型での階乗計算
【Java】1000000!のような巨大な階乗計算をする [任意精度整数(多倍長整数)]
【Java】任意精度整数(多倍長整数)で演算する

■参考になる書籍

関連記事

category: Java

thread: プログラミング

janre: コンピュータ

tag: 算術関数動的計画法最短経路練習問題

tb: 0 cm: --

【Java】ナップサック問題(knapsack)[動的計画法]

2015/04/28 Tue [edit]

　最強最速アルゴリズマー養成講座という本の P.177～186 にナップサック問題(knapsack)の解法が載っているのだが(C#,C++,Java 版共に)、品物データの並びを変えると答えが間違って出てくることがあるので、分析してみたらいくつかバグが見つかった。P.180～181 のメモ化再帰でも「return -1;」となっている部分は「return -ps[n - 1];」にしないとバグると思うが、P.184～185 のボトムアップの動的計画法の方は「入れない」処理がされてないので、やはり並びによってはバグる。

　正誤表には載ってないね。検索しても修正版など見つからなかったので、結局 Wikipedia を見ながら書き直したら上手くいった。

●P.184～185 ナップサック問題の動的計画法での解法 O(nW)

int[] ws = {3, 4, 1, 2, 3};  //重さ
int[] ps = {2, 3, 2, 3, 6};  //価値
int maxw = 10;  //重量制限
int[][] dp = new int[ws.length + 1][maxw + 1];
int ret = 0;

for (int i = 0; i < ws.length; i++) {
    for (int j = 0; j <= maxw; j++) {
        if (ws[i] <= j) {
            dp[i + 1][j] = Math.max(dp[i][j], dp[i][j - ws[i]] + ps[i]);
            ret = Math.max(ret, dp[i + 1][j]);
        } else {
            dp[i + 1][j] = dp[i][j];
        }
    }
}
//※答えは dp[ws.length][maxw] でも良い

　C#, C++ 版も同じようなものなので、書き換えれば上手くいくと思う。

　ちなみに問題の内容としては以下の感じになる。

それぞれの重さが ws、価値が ps の品物があります。これらから好きなものを選んで、ナップサックに詰め込みたいと思います。ただしこのナップサックは全部で 10kg までしか入りません。ナップサックに詰め込める品物の価値の和の最大値を求めなさい。

　たぶん元のコードに品物を「入れない」処理

dp[i + 1][j] = Math.max(dp[i + 1][j], dp[i][j]);

が抜けているので、それ以前の品物の引き継ぎ状態が消えているのだと思う。配列を出力してみればいくつか抜けているね。あとは不等号やループ変数を ps[j] → ps[i] に修正すれば、上手く動くね。以下のような感じかな。

int[] ws = {3, 4, 1, 2, 3};
int[] ps = {2, 3, 2, 3, 6};
int maxw = 10;
int[][] dp = new int[ws.length + 1][maxw + 1];
int ret = 0;

for (int i = 0; i < ws.length; i++) {
    for (int j = 0; j <= maxw; j++) {
        dp[i + 1][j] = Math.max(dp[i + 1][j], dp[i][j]);
        if (j + ws[i] <= maxw) {
            dp[i + 1][j + ws[i]] = Math.max(dp[i + 1][j + ws[i]], dp[i][j] + ps[i]);
            ret = Math.max(ret, dp[i + 1][j + ws[i]]);
        }
    }
}
//※答えは dp[ws.length][maxw] でも良い

　アルゴリズムの解説は、グラフや配列の中身を詳しく書いてあるブログなどがあるのでそちらに譲ろう。私流の解釈としては「それぞれ 1kg, 2kg, 3kg,..., 10kg まで入る器を用意して、次々と品物を『入れるとき』と『入れないとき』(現状のまま)を比べて、価値が大きくなる方だけを残すことを繰り返すと、最後には最大の価値がわかる」と考えればわかりやすい気がする。

　そう考えれば、経過を保存している２次元配列 dp[][] の代わりに、１次元配列を２つ用意して、コピー元とコピー先をひっくり返す手もあるね。コードを見ても直前のインデクスしか参照してないし、２列だけあれば事足りる。データが大量にあるとき、メモリの節約にはなるかもね。ちなみに最大値を保存する ret は、最後の結果の配列の一番右になるので必要ない。また重さ０の列は最後までずっと０なので、ループも１からで良いかも知れない（※Java の場合、配列は０で初期化されるので）。

●経過保存の配列を１次元×２つにして、ちょっぴりメモリ節約版？

int[] src = new int[maxw + 1];    //dp[][] の代わり
int[] dst = new int[maxw + 1];    //dp[][] の代わり
int[] tmp;    //スワップ用

for (int i = 0; i < ws.length; i++) {
    for (int j = 1; j <= maxw; j++) {
        if (j < ws[i]) {
            dst[j] = src[j];
        } else {
            dst[j] = Math.max(src[j], src[j - ws[i]] + ps[i]);
        }
    }
    //参照をスワップ(C言語系ならポインタをスワップ)
    tmp = src;
    src = dst;
    dst = tmp;
}
//最後の結果は src[maxw] になる(スワップしているため)