ヽ|∵|ゝ(Fantom) の 開発blog? ホーム »練習問題
このページの記事一覧

【Java】拡張ユークリッドの互除法  


 「ユークリッドの互除法」は「自然数 a, b の最大公約数」を求めるものだったが、「拡張ユークリッド互除法」は「ax + by = gcd(a, b) となる a, b の最大公約数を求める」ものになる。簡単な方程式の整数解を求めるのにも使える。


 ググッてみると C言語で書かれている関数で引数が extgcd(a, b, &x, &y) ['&'は参照渡しとなる] となっているものが多いが、Java はプリミティブ型で参照渡しができないので、ループ版(非再帰版)の方は英語版 Wikipedia に載っているように、引数を extgcd(a, b) にして作ってみよう。無理矢理参照渡しのようにするために、要素が1つの配列を用いる方法もあるが、それは再帰版でのサンプルを見て欲しい。ループ版も同じように書き換えることもできるが、3つ以上の係数に対応したいときなどには不便なので、用途によって使い分けても良いかも知れない。

(参考)
ユークリッドの互除法/拡張ユークリッドの互除法
拡張版ユークリッドの互除法
拡張ユークリッドの互除法 (Extended Euclidean Algorithm)
Extended Euclidean algorithm

●拡張ユークリッドの互除法(非再帰版)
/**<h1>拡張ユークリッド互除法</h1>
* <p>ax + by = gcd(a, b) となる a, b の最大公約数と解 x, y を求める。</p>
* @param a : 数値1(>0)
* @param b : 数値2(>0)
* @return<b>int[]</b> : {[0]:gcd, [1]:x, [2]:y}:最大公約数(なし=1 [互いに素])と解 x, y
*/
public static final int[] extgcd(int a, int b) {
int x0 = 1, x1 = 0;
int y0 = 0, y1 = 1;

while (b != 0) {
int q = a / b;
int r = a % b;
int x2 = x0 - q * x1;
int y2 = y0 - q * y1;

a = b; b = r;
x0 = x1; x1 = x2;
y0 = y1; y1 = y2;
}

return new int[]{a, x0, y0};
}


//メインでは...
int a = 12;
int b = 21;
int[] res = extgcd(a, b);
int g = res[0];
int x = res[1];
int y = res[2];
System.out.println("g = " + g + ", x = " + x + ", y = " + y);
System.out.println(a + " * " + x + " + " + b + " * " + y + " = " + g);
int sum = a * x + b * y;
System.out.println(sum == g);

g = 3, x = 2, y = -1
12 * 2 + 21 * -1 = 3
true

 出力の意味は a = 12, b = 21 のときに、12 * 2 + 21 * -1 = 3 [ax + by = gcd(a, b)] となる解 x = 2, y = -1 が存在することを表す。実はこのような解は常に存在するので、このことを利用して解く問題などもある。

 また、最大公約数がないときは 1 となるが、このとき「互いに素」であることがわかる(お互いの共通の因子を持たない)。しばしばこの性質を利用する解法も見かけるので、覚えておいた方が良いだろう。

 もう1つ、再帰でも書いてみよう。冒頭に書いた「無理矢理参照渡しのようにするために要素が1つの配列を用いる方法」の例でもある。表記が複雑になるので、少し不便だが、C言語などからの移植する際の手法の1つとして覚えておくと役に立つこともあるだろう。

●拡張ユークリッドの互除法(再帰版)
/**<h1>拡張ユークリッド互除法</h1>
* <p>ax + by = gcd(a, b) となる a, b の最大公約数と解 x, y を求める。</p>
* @param a : 数値1(>0)
* @param b : 数値2(>0)
* @param x : 解 x (参照渡し用 new int[1] で作る)
* @param y : 解 y (参照渡し用 new int[1] で作る)
* @return<b>int</b> : 最大公約数(なし=1 [互いに素])
*/
public static final int extgcd(final int a, final int b, final int[] x, final int[] y) {
int g = a;
if (b != 0) {
g = extgcd(b, a % b, y, x);
y[0] -= (a / b) * x[0];
} else {
x[0] = 1;
y[0] = 0;
}
return g;
}


//メインでは...
int a = 12;
int b = 21;
int[] x = new int[1]; //参照渡しにするため、要素1個で作る
int[] y = new int[1]; //参照渡しにするため、要素1個で作る
int g = extgcd(a, b, x, y);
System.out.println("g = " + g + ", x = " + x[0] + ", y = " + y[0]);
System.out.println(a + " * " + x[0] + " + " + b + " * " + y[0] + " = " + g);
int sum = a * x[0] + b * y[0];
System.out.println(sum == g);

g = 3, x = 2, y = -1
12 * 2 + 21 * -1 = 3
true

 ループ版と再帰版とで実行速度を比べてみると、差はほとんどないので、どちらを使っても良いだろう。


 せっかくなので、上記の extgcd() を使って練習問題を解いてみよう。

●双六(すごろく) (参考) アリ本 P.108 [解答例は載ってない]
双六(すごろく)があります。マス目には整数が書かれています。マス目は前後に無限に続いており、
0のマス目がスタートで、1のマス目がゴールなのですが、サイコロの目には a, b, -a, -b の4つの整数しかありません。
そのため、a, b の値によってはゴールに辿り着けない場合があります。

 ┌─┬─┬─┬─┬─┬─┬─┬─┬─┐
…│-4│-3│-2│-1│ 0│ 1│ 2│ 3│ 4│…
 └─┴─┴─┴─┴─┴─┴─┴─┴─┘

4つの目をそれぞれ何回出せばゴールに辿り着けますか?
複数解が存在する場合は、どれか1つを出力しなさい。解がなければ -1 を出力しなさい。

制約
・1 <= a, b <= 109

入力
a = 4, b = 11

●拡張ユークリッド互除法での解法
int a = 4;    //入力
int b = 11; //入力

int[] res = extgcd(a, b);
int g = res[0];
int x = res[1];
int y = res[2];

if (g != 1) {
System.out.println(-1);
} else {
String ans = "" + ((x >= 0) ? x : 0);
ans += " " + ((y >= 0) ? y : 0);
ans += " " + ((x < 0) ? -x : 0);
ans += " " + ((y < 0) ? -y : 0);
System.out.println(ans); //答え
}

3 0 0 1

 答えは左から a, b, -a, -b に対応させて、右へ 4(a)×3回、左へ 11(b)×1回 進めば、1 (gcd(a, b)) のマス目(ゴール)に辿り着けることを表している。ちなみに a = 11, b = 4 と入れ替えて試してみると、答えは「0 3 1 0」となり意味は同じとなる。では逆に辿り着けない例はどうなるだろうか。例えば gcd が2となる a, b を与えれば(両方偶数など)、1には辿り着けない。全探索で見つけるまでもなく一瞬で求められるのは非常に有り難い(笑)。


(関連記事)
【Java】3つ以上の最大公約数を求める(ユークリッドの互除法②)
【Java】最大公約数・最小公倍数を求める(ユークリッドの互除法)
【Java】約数を求める
【Java】素因数分解をする


■参考になる書籍


スポンサーサイト

category: Java

thread: プログラミング

janre: コンピュータ

tag: 算術関数  約数  練習問題 
tb: 0   cm: --

【Java】最大公約数・最小公倍数を求める(ユークリッドの互除法)  


 以前ちらりと書いたが、「エラトステネスの篩」と同じように、コンピュータが存在する以前からある「最大公約数」(Greatest Common Divisor) を求めるアルゴリズムに「ユークリッドの互除法」というものがある。最古のアルゴリズムとも言われているが、コンピュータ時代になっても健在で、簡単でしかも高速なアルゴリズムとしてもよく利用されている。


 概要は「正の整数 a, b (a >= b) について「a の b による剰余を r とすると、 a と b との最大公約数は b と r との最大公約数に等しいという性質が成り立つ。この性質を利用して、 b を r で割った剰余、 除数 r をその剰余で割った剰余、と剰余を求める計算を逐次繰り返すと、剰余が 0 になった時の除数が a と b との最大公約数となる」(Wikipedia) となっているが、Wikipedia の右端にあるアニメーションを見た方がわかりやすい気がする(笑)。

 書き方もループで書いたもの再帰で書いたものをよく見かけるが、比較のためにも両方書いてみよう。

●最大公約数を求める [ユークリッドの互除法](非再帰版)
/**<h1>最大公約数を求める(Greatest Common Divisor)[ユークリッドの互除法]</h1>
* @param a : 数値1(>0)
* @param b : 数値2(>0)
* @return<b>int</b> : 最大公約数 (なし=1 [互いに素])
*/
public static final int gcd(int a, int b) {
//a > b にする(スワップ)
if (a < b) {
int tmp = a;
a = b;
b = tmp;
}

//ユークリッドの互除法
int r = -1;
while (r != 0) {
r = a % b;
a = b;
b = r;
}

return a; //b には r=0 の値が入るため、a を返す
}


//メインでは...
int a = 1920;
int b = 1080;
int d = gcd(a, b);
System.out.println("d = " + d); //最大公約数
System.out.println(a + "/" + b + " = " + (a / d) + "/" + (b / d)); //約分に利用

d = 120
1920/1080 = 16/9

 例では 1920 と 1080 の最大公約数 d は 120 であり、それぞれを最大公約数(120)で割ると、約分と同じ意味になることを表している。1920x1080 はアスペクト比で有名なので、16 : 9 (16/9) が正しいのは検算しなくてもわかるだろう。

 関数冒頭で、a, b の大小関係を調べてスワップしているが、ユークリッドの互除法は「a, b が正の値で a >= b とするとき」となっているので、既知ならば削除しても構わない。

 また、最大公約数がないときは 1 となるが、このとき「互いに素」であることがわかる(お互いの共通の因子を持たない)。実はしばしばこの性質を利用する解法も見かけるので、覚えておいた方が良いだろう。

 もう1つ、ユークリッドの互除法は再帰の例としても挙げられるが、確かに非常に簡潔に書けるので、あまり負荷のかからない処理なら、こちらでも構わないだろう。

●最大公約数を求める [ユークリッドの互除法](再帰版)
/**<h1>最大公約数を求める(Greatest Common Divisor)[ユークリッドの互除法](再帰版)</h1>
* @param a : 数値1(>0)
* @param b : 数値2(>0)
* @return<b>int</b> : 最大公約数 (なし=1 [互いに素])
*/
public static final int gcd(final int a, final int b) {
return (b == 0) ? a : gcd(b, a % b);
}


//メインでは...
int a = 1920;
int b = 1080;
int d = (a >= b) ? gcd(a, b) : gcd(b, a); //a >= b にする
System.out.println("d = " + d); //最大公約数
System.out.println(a + "/" + b + " = " + (a / d) + "/" + (b / d)); //約分に利用

d = 120
1920/1080 = 16/9

 こちらは随分簡潔になっているが、非再帰版は a, b の大小関係はチェックしてないので、注意しよう。心配ならサンプルのようにメインコードで調べれば良い。

 また、ループ版と再帰版とで実行速度を比べてみると、やはりループ版の方がわずかに速いが、それほど差はないので、計算量の小さいものなら問題なく使えるだろう。


 上の例では gcd() を約分に利用しているが、他にも色々と応用できる。例えば「線分上の格子点の個数」を求めてみよう。

●線分上の格子点の個数 (参考) アリ本 P.107 [解答例は載ってない]
平面の2つの格子点 P1 = (x1, y1), P2 = (x2, y2) が与えられます。
線分 P1, P2 上には P1, P2 以外にいくつの格子点が存在しますか?
(格子点とは) 座標平面(or 座標空間)において,各成分が全て整数であるような点。

制約
・-109 <= x1, x2, y1, y2 <= 109

入力例
P1 = (1, 11), P2 = (5, 3)

● 線分上の格子点の個数を求める gcd() での解法
int x1 = 1, y1 = 11;    //P1
int x2 = 5, y2 = 3; //P2

int a = Math.abs(x1 - x2);
int b = Math.abs(y1 - y2);
int g = 1;

if (x1 == x2 && y1 == y2) {
//0
} else if (x1 == x2) {
g = b;
} else if (y1 == y2) {
g = a;
} else {
g = gcd(a, b);
}
System.out.println(g - 1); //答え

//以下は具体的な座標の出力
int x = x1, y = y1, dx = 0, dy = 0;
if (x1 <= x2) {
dx = (x2 - x1) / g;
dy = (y2 - y1) / g;
} else {
x = x2; y = y2;
dx = (x1 - x2) / g;
dy = (y1 - y2) / g;
}
for (int i = 1; i < g; i++) {
if (i > 1) {
System.out.print(", ");
}
System.out.print("(" + (x + i * dx) + ", " + (y + i * dy) + ")");
}
System.out.println();

3
(2, 9), (3, 7), (4, 5)

 単純に (x1, y1)-(x2, y2) の線分上でかつ整数となる点を全て調べても良いが、線分が大きくなると計算量が間に合わなくなる。その点、gcd() なら非常に軽い。gcd() は0で除算できない(P1, P2 が同じ点、または x 軸, y 軸のどちらかが同じ)ことだけ注意しよう。

- - - - - - - - - - - - - - - - - - 
●最小公倍数を求める

 ついでにもう1つ、最小公倍数(Least Common Multiple)も求めてみよう。

 Wikipedia にも書いてあるように「正の整数a, bに対して、最大公約数gcd(a, b)と最小公倍数lcm(a, b)との間には、gcd(a, b)×lcm(a, b) = ab という関係がある」となっているので、「gcd()」ができていれば、非常に簡単に作れる。

●最小公倍数を求める
/**<h1>最小公倍数を求める(Least Common Multiple)</h1>
* @param a : 数値1(>0)
* @param b : 数値2(>0)
* @return<b>int</b> : 最小公倍数
*/
public static final int lcm(final int a, final int b) {
return a * b / gcd(a, b);
}


//メインでは...
int h = 16;
int v = 9;
int m = lcm(h, v);
System.out.println("m = " + m);

m = 144

 こちらも「gcd()」のおかげで随分と簡潔だ。内容は式を変形しただけのものに過ぎないので、説明はいらないだろう。

 「ユークリッド互除法」には「ax + by = gcd(a, b) となる a, b の最大公約数を求める」という「拡張ユークリッド互除法」というものもある。解となる整数 x, y の組を見つけることができるものなので、数学的な問題には役に立つことがあるかも知れない。実装は gcd() とそれほど変わらないので。覚えておいても損はないだろう。


■約分・最大公約数・最小公倍数 計算機
ユークリッドの互除法で約分・最大公約数(gcd)・最小公倍数(lcm)を求める(JavaScript版)
元の分数を入力:
/



約分された分数:
/

最大公約数(gcd):
最小公倍数(lcm):


(関連記事)
【Java】3つ以上の最大公約数を求める(ユークリッドの互除法②)
【Java】拡張ユークリッドの互除法
【Java】約数を求める
【Java】素因数分解をする
【Java】素数判定①(試し割り法)


■参考になる書籍


category: Java

thread: プログラミング

janre: コンピュータ

tag: 算術関数  約数  練習問題  アルゴリズム 
tb: 0   cm: --

【Java】式を逆ポーランド記法(後置記法)で計算する  


 今度は前回「逆ポーランド記法(後置記法)への変換」で変換された文字列を実際に計算して数値にする関数を作ろう。

 といってもコードは非常にシンプルなので驚くかも知れない。これが逆ポーランド記法がコンピュータと相性が良いと言われる所以かも知れない。ただひたすらスタックからオペランドと演算子を取り出し、その結果を再びスタックに積むという動作を繰り返すだけだ。

●逆ポーランド記法(Reverse Polish Notation)[=後置記法:postfix notation]を計算する
import java.util.ArrayDeque;
import java.util.Deque;

/**<h1>逆ポーランド記法(Reverse Polish Notation)[=後置記法:postfix notation]を計算する</h1>
* <p>スペースで区切った "10 20 30 * +" のような形式。連続した空白は1つ分の空白として処理される。
* <br>四則演算(+, -, *, /)と括弧のみ。</p>
* @param rpn : 逆ポーランド記法での計算式
* @return<b>int</b> : 計算された値
*/
public static final int parseRPN(final String rpn) {
final Deque<Integer> stack = new ArrayDeque<Integer>();

String s = rpn.replaceAll("\\s+", " "); //連続した空白文字をスペース1つ分に置き換える
final String[] arr = s.trim().split(" ");

for (String e : arr) {
if ('0' <= e.charAt(0) && e.charAt(0) <= '9') { //数字
stack.push(Integer.parseInt(e));
} else {
int a = stack.pop();
int b = stack.pop();
if (e.equals("*")) {
stack.push(b * a);
} else if (e.equals("/")) {
stack.push(b / a); //div/0 に注意
} else if (e.equals("+")) {
stack.push(b + a);
} else if (e.equals("-")) {
stack.push(b - a);
}
}
}

return stack.pop();
}


//メインでは...
String rpn = "10 20 + 30 40 - * 50 2 / +";
System.out.println(parseRPN(rpn));

-275

 注意点は前回と同じなので「式を逆ポーランド記法(後置記法)変換する」の記事を参照して欲しい。

 計算される値はとりあえず int 型になっているので、long で算出したいなら、型を修正するのも良いだろう。ゼロ除算(エラー)や剰余などは付けてないので、必要あれば追加しても良いだろう。特に負の除算や剰余に関しては、言語によって仕様が違うので、実装前に予めパースする式の符号などを確認しておく必要がある。以下に参考ページも載せておこう。参考は「負の剰余」がテーマになっているが、商と剰余は2つで1セットのようなものなので、どちらも確認しておいた方が良いだろう。

(参考)
負数の剰余を計算してはならない
こんなプログラムはいやだ: 負の剰余
負の剰余

 また、前回の「逆ポーランド記法への変換」と今回の「逆ポーランド記法を計算」の両方を同時に行うと、簡単な数式をパースして計算する関数にもなる。これは簡易的な eval() 関数やインタプリタなどにも使える。実は元々はそのために作ったものだったりする(笑)。コードを見てもらえばわかると思うが、本当にそのまま2つの関数を混ぜあわせただけだ。

●数式をパースして計算する
import java.util.ArrayDeque;
import java.util.Deque;
import java.util.HashMap;
import java.util.Map;

//※rpnRank は前回の Map をそのまま利用なので割愛。

/**<h1>数式をパースして計算する</h1>
* <p>逆ポーランド法で計算する。使える演算子は四則演算と括弧のみ。
* <br>スペースを含まない "(10+20)*30-(40+50)/2" のような式。空白文字は全て除去される。</p>
* @param expression : 数式の文字列
* @return<b>int</b> : 計算値
*/
public static final int parseExpression(final String expression) {
final Deque<Character> stack = new ArrayDeque<Character>();
final Deque<Integer> val = new ArrayDeque<Integer>();

String s = expression.replaceAll("\\s+", ""); //空白文字を取り除く
s = "(" + s + ")"; //末尾に")"を付けることで、最後にスタックを吐き出させる
final int len = s.length();

String tmp = "";
for (int i = 0; i < len; i++) {
char c = s.charAt(i);
if ('0' <= c && c <= '9') {
tmp += c; //数字1文字ずつのため
} else {
if (!tmp.equals("")) {
val.push(Integer.parseInt(tmp));
tmp = "";
}

while (!stack.isEmpty() && rpnRank.get(stack.peek()) >= rpnRank.get(c) && stack.peek() != '(') {
char e = stack.pop();
int a = val.pop();
int b = val.pop();
if (e == '*') {
val.push(b * a);
} else if (e == '/') {
val.push(b / a); //div/0 に注意
} else if (e == '+') {
val.push(b + a);
} else if (e == '-') {
val.push(b - a);
}
}
if (c == ')') {
stack.pop(); //'('
} else {
stack.push(c);
}
}
}

return val.pop();
}


//メインでは...
String expr = "(10+20)*(30-40)+50/2";
System.out.println(parseExpression(expr));

-275

 解説は同じものなので特に必要はないだろう。

 せっかくなので、この関数を使って練習問題を解いてみよう。問題は「プログラマ脳を鍛える数学パズル」から抜粋させて頂いた。この本では主に Ruby, JavaScript, C で解答ソースが書かれているが、アルゴリズムの本質というものは言語とは特に関係ないので(考え方の方法なので)、理解してしまえば色々応用できる。今回は Ruby ソースを Java に移植したものと考えても良いだろう。

●Q42 1つの数字で作る 1234 (IQ:100 目標時間:30分) ★★★
1つの数字だけで、四則演算によってある数を作ることを考えます。
例えば「1000」を作るとき、「1」だけを7個使って「1111-111」で表現できます。
「8」だけなら8個使って「8+8+8+88+888」、「9」だけであれば5個使って「9÷9÷999」で1000になります。
使用できるのは加減乗除だけとし、括弧は使用不可とします。また「÷」は整数の商のみとします(端数切捨て)。

「1234」を1つの数字だけで、できるだけ少ない個数で表現するとき、最も少ない個数で表現できる数字を求め、
その式をすべて答えて下さい。

●上記の parseExpression() を使った Java での解法
import java.util.ArrayDeque;
import java.util.Deque;
import java.util.HashMap;
import java.util.Map;

public class Q42 {
public static void main(String[] args) throws Exception {
new Q42();
}

public Q42() throws Exception {
// n = 1000; //9/9+999, 999+9/9
n = 1234; //99999/9/9
// n = 4321; //6/6+66*66-6*6, 6/6-6*6+66*66, 66*66+6/6-6*6, 66*66-6*6+6/6
// n = 2016; //9+9+999+999, 9+999+9+999, 9+999+999+9, 999+9+9+999, 999+9+999+9, 999+999+9+9

solve();
}

static final String[] op = {"+", "-", "*", "/", ""};
int n; //作る数
boolean found = false;

void solve() {
int len = 1; //使用する桁数
while (!found) {
for (int i = 1; i <= 9; i++) {
check(len, String.valueOf(i), i);
}
len++;
}
}

//remain:桁残り, expr:ここまでの式, num:使用する数字
void check(int remain, String expr, int num) {
if (remain == 0) {
if (parseExpression(expr) == n) { //※上記の関数を使う
System.out.println(expr);
found = true;
}
} else {
for (String o : op) { //演算子 or ""
check(remain - 1, expr + o + String.valueOf(num), num);
}
}
}

//※parseExpression() は割愛(上記のものを使う)
}

99999/9/9

 paiza.io(オンラインIDE)で試してみると、問題の「1234」の場合、約 500ms くらいなので、ギリギリといったところか(笑)。「4321」や「2016」のように答えが多いものは 1s を超えるので、オンラインジャッジでは少しマズイ(1つだけで良いのなら、探索を中断すれば問題ない)。

 ちなみに「ScriptEngine」で外部言語として eval() を使うには以下のようなコードになる。

●「ScriptEngine」で JavaScript の eval() を使って数式を計算する
import javax.script.ScriptEngine;
import javax.script.ScriptEngineManager;

public class Main {
public static void main(String[] args) throws Exception {
new Main();
}

public Main() throws Exception {
ScriptEngineManager manager = new ScriptEngineManager();
ScriptEngine engine = manager.getEngineByName("JavaScript");

String expr = "(10+20)*(30-40)+50/2";
System.out.println(engine.eval(expr));;
}
}

-275

 ただし、非常に重いので(数式2つ以上で 1s 超える(paiza.io))全探索で利用するのは難しいだろう。

 また、この「parseExpression()」を使えば、「Q02 数列の四則演算」の問題も同じような方法で解けるので挑戦してみるのも良いだろう。

●Q02 数列の四則演算 (IQ:70 目標時間:10分) ★
並んでいる数字の各桁の間に四則演算の演算子を入れて計算した結果「元の数を桁を逆から並べた数字」
と同じになるものを考えます(演算子は入れない場所があっても構いませんが、最低でも1つは入れるものとします)。
例えば、100~999 (3桁の数字) の場合、以下の 3つがあります。

351 → 3×51 = 153
621 → 6×21 = 126
886 → 8×86 = 688

1000~9999 (4桁の数字) のうち、上記の条件を満たす数を求めて下さい。

[答え] 5931 (→ 5×9×31 = 1395)
[>>クリックでコードを見る]


※なお、問題文などは以下の書籍からお借りしている(ただし、書籍の解答例は Ruby, JavaScript, C などで書かれているので購入の際はご注意)。パズルっぽい問題の解法パターンを勉強するには良いかも。




(関連記事)
【Java】式を逆ポーランド記法(後置記法)に変換する


category: Java

thread: プログラミング

janre: コンピュータ

tag: 算術関数  練習問題  アルゴリズム 
tb: 0   cm: --

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


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



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

1.頂点(一番上)に「1」を書き、ここをスタート地点とする。

2.スタート地点から左右(矢印の方向:斜め下)に向かって進んでいく。その際、交差する地点はその1つ前の値を合計する。これを繰り返す。

 これは道に例えると「網の目のように必ず2つの分岐点がある道を進むとき、その地点に最短で辿り着く方法は何通りあるか?」という数になる。つまり、地点Aに辿り着く方法が2通り、地点Bに辿り着く方法が1通りあるとき、その地点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(n2)]
/**<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)]と比べたメリットは、計算が単純な加算しかしていないため、オーバフローが少しばかりしにくい所だろうか。また、大量に計算するなら、1つパスカルの三角形を作っておいて[O(n2)]、配列から値を取り出すだけにすれば[O(1)]、実行負荷はほとんど無くなる(例えば、30個の組合せを1000回計算するなら、漸化式では 30[O(n)]x1000回 = 30000回演算、パスカルの三角形なら 30x30[O(n2)]=900 +1000[O(1)]回 = 1900回演算で、計算量が少なくなるので速い)。しかし単体で計算するなら、毎回配列を作る[O(n2)]ので、実行速度は遅くなるので注意しよう。


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

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

//パスカルの三角形での解法
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; //グー、チョキ、パーの3手でそれぞれ勝てるため、3倍する
double d = 1.0 - sum / Math.pow(3, n); //あいこ = 全事象 - 勝敗がつく事象 [=余事象]
System.out.println(d); //0.4814814814814815 (13/27)

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

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

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


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

//あいこになる確率を一般化したものでの解法
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人でじゃんけんを1回するとき、あいこになる確率 = 0.33333333333333337
3人でじゃんけんを1回するとき、あいこになる確率 = 0.33333333333333337
4人でじゃんけんを1回するとき、あいこになる確率 = 0.4814814814814815
5人でじゃんけんを1回するとき、あいこになる確率 = 0.6296296296296297
6人でじゃんけんを1回するとき、あいこになる確率 = 0.7448559670781894
7人でじゃんけんを1回するとき、あいこになる確率 = 0.8271604938271605
8人でじゃんけんを1回するとき、あいこになる確率 = 0.8838591678097851
9人でじゃんけんを1回するとき、あいこになる確率 = 0.922267946959305
10人でじゃんけんを1回するとき、あいこになる確率 = 0.9480770207793527
11人でじゃんけんを1回するとき、あいこになる確率 = 0.9653508103439516
12人でじゃんけんを1回するとき、あいこになる確率 = 0.9768892501707621
13人でじゃんけんを1回するとき、あいこになる確率 = 0.9845890700943284
14人でじゃんけんを1回するとき、あいこになる確率 = 0.9897247922786035
15人でじゃんけんを1回するとき、あいこになる確率 = 0.9931494433687528
16人でじゃんけんを1回するとき、あいこになる確率 = 0.9954328228623964
17人でじゃんけんを1回するとき、あいこになる確率 = 0.9969551687804513
18人でじゃんけんを1回するとき、あいこになる確率 = 0.9979700970332521
19人でじゃんけんを1回するとき、あいこになる確率 = 0.9986467261931519
20人でじゃんけんを1回するとき、あいこになる確率 = 0.9990978157413181

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


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

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

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


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

 実はこれも高校数学の教科書に載っているのだが、解法を簡単に書いてみると、下へ進む道を a 、右に進む道を b として1列に並べると、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(2hw) くらいの計算量はかかってしまうので、非常に重くなる。あと、マスが 100 x 100 = 10000 を超えると Stack Overflow を起こすこともあるので気をつけよう。このパスカルの三角形を改造した解法(動的計画法)なら、計算量は O(hw) でしかない。

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


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


■参考になる書籍


category: Java

thread: プログラミング

janre: コンピュータ

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

【Java】組合せの数を求める  


 「階乗」「順列」とやったので、「組合せ」もやってみよう。

 「組合せ」とは簡単に言えば「n 個ものもから、r 個取り出す方法で、並べ方を問題にしない数」である。「8人の中から3人選び出してチームを作る方法は何通りあるか?」のような数を計算するのに使う。並べ方を考慮する場合は「順列」になる。

 組合せの計算は以下の式のようになるが、これも階乗計算同様、桁あふれのオーバフローを起こす可能性が高いので、一番下の「漸化式」を使おう。こちらの式を使うと多少オーバフローを防ぐことができる。また、ちょっとした計算として使うなら、オーバーフローの心配のない Ruby 版を使う手もある。開発環境がなくてもオンラインIDEで十分だろう。
 


●組合せの計算
/**<h1>n 個のものから r 個取り出した組合せの数を求める</h1>
* <p>漸化式 : nCr = (n - r + 1) / r * nCr-1 を利用する。</p>
* @param n : すべての数 (1~61) (負の値はエラー)
* @param r : 取り出す数 (1~30) (負の値はエラー)
* @return<b>long</b> : 組合せの数
*/
public static final long combination(final int n, int r) {
if (n < 0 || r < 0 || n < r) {
throw new ArithmeticException("n = " + n + ", r = " + r);
}

r = Math.min(r, n - r); //nCr = nCn-r
if (r == 1) {
return n;
}

long sum = 1;
for (int i = 1; i <= r; i++) {
//漸化式 : nCr = (n - r + 1) / r * nCr-1 を利用
sum = sum * (n - i + 1) / i; //漸化式(r(i)で割り切るのに、先に nCr-1 が必要なため、 sum *= ~ としない)
}
return sum;
}


//メインでは...
int n = 20;
int r = 10;
long l = combination(n, r);
System.out.println(l);

 引数は2つとり、負の値では色々都合悪いので例外処理も入れた。オーバフローは n = 61 のとき r = 30 までは可能だが、心配ならループ内に負の値のチェックを入れるのも良いだろう(もちろんその分、実行速度は落ちる)。結果を double にすれば、指数形式にはなるが、上から16桁までは計算できる(17桁以降は丸められ、0 で埋まる)

 組合せの nCr = nCn-r の性質から、r と n-r の小さい方をとったり、r = 1 のときは計算する必要ないので、条件を1つ付け加えたりしているが、少しでも計算を速くするためなので、コードを短くしたいなら取り除いても構わない(結果は変わらない)。

 ちなみに、「n 種類のものから、重複して取り出すことを許して、r 個の集まりを作る」という「重複組合せ」というものもあるが、式を比較してみると、n の値に (n + r - 1) を代入したものと同じになっていることがわかる。だから上記の関数を使って、combination(n + r - 1, r) で求められることもわかる。


 もう1つ、オーバフロー対策として任意精度整数(多倍長整数)を使って書いてみよう。実行速度はプリミティブ型と比べて約35倍ほど劣るが、桁あふれを防げるので検算には役に立つだろう。また、再帰で書く方法はただでさえ重い計算なので、避けたほうが良いだろう(約2倍の時間がかかる)[→計算時間は階乗計算の記事を参考に]。
 
●任意精度整数での組合せの計算
import java.math.BigInteger;

public static final BigInteger combination(final int n, int r) {
if (n < 0 || r < 0 || n < r) {
throw new ArithmeticException("n = " + n + ", r = " + r);
}

r = Math.min(r, n - r); //nCr = nCn-r
if (r == 1) {
return new BigInteger(String.valueOf(n));
}

BigInteger sum = BigInteger.ONE;
for (int i = 1; i <= r; i++) {
//漸化式 : nCr = (n - r + 1) / r * nCr-1 を利用
sum = sum.multiply(new BigInteger(String.valueOf((n - i + 1)))).divide(new BigInteger(String.valueOf(i)));
}
return sum;
}


//メインでは...
int n = 60;
int r = 30;
BigInteger b = combination(n, r);
System.out.println(b);


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

4枚のコインを投げた時、そのうち2枚が表で2枚が裏になる確率はどれくらいか?

//組み合わせでの解法
double d = (double)combination(4, 2) / Math.pow(2, 4);
System.out.println(d); //0.375 (= 3/8)

 答えは 37.5% (= 3/8) である。「えっ?50% (1/2) じゃないの?」と思った人は気をつけよう。それはコイン4枚のうち2枚という数字で 2/4 = 1/2 のような思い込みをしているだけだ。

 簡単に解法を書いてみると、コインは表裏の2通りで4枚あるから、すべての可能性は 2 x 2 x 2 x 2 = 16 通りである。うち2枚が表(残りは必然的に裏)というのは組合せの 4C2 = 6 通りだから、確率は (起こる事象 / 全事象) で 6/16 = 3/8 = 0.375 (37.5%) となる。

 2進数がわかる人なら(プログラマーなら知ってると思うが)、表を1、裏を0としたとき、その並びは 0000 ~ 1111 つまり 0 ~ 15 (16 個) というのはすぐわかるし、そのうち2つのビットが1になっているものは、0011, 0101, 0110, 1001, 1010, 1100 の6つだから、確率は 6/16 = 3/8 となるのもわかるだろう。

 この問題なら数が少ないので、自分で数えてもすぐにわかると思うが、では「50枚のコインを投げて、25枚表、25枚裏になる確率は?」と問われたらどうだろうか?もちろん 1/2 ではない(笑)。とたんに手計算ではやってられなくなるだろう。その場合、先ほどのコードを変数に置き換えてみると、

int n = 50;
int r = 25;
double d = (double)combination(n, r) / Math.pow(2, n);
System.out.println(d); //0.11227517265921705

のようになる(答え 11.2%)。

 「あれ?4枚のときより少なくなっているな…」と気づいた人はスルドイ。試しに 4枚~50枚 でその半分の枚数が表になる確率を計算してみよう。

4枚のうち 2 枚が表の確率 = 0.375
6枚のうち 3 枚が表の確率 = 0.3125
8枚のうち 4 枚が表の確率 = 0.2734375
10枚のうち 5 枚が表の確率 = 0.24609375
12枚のうち 6 枚が表の確率 = 0.2255859375
14枚のうち 7 枚が表の確率 = 0.20947265625
16枚のうち 8 枚が表の確率 = 0.196380615234375
18枚のうち 9 枚が表の確率 = 0.1854705810546875
20枚のうち 10 枚が表の確率 = 0.17619705200195312
22枚のうち 11 枚が表の確率 = 0.16818809509277344
24枚のうち 12 枚が表の確率 = 0.1611802577972412
26枚のうち 13 枚が表の確率 = 0.15498101711273193
28枚のうち 14 枚が表の確率 = 0.14944598078727722
30枚のうち 15 枚が表の確率 = 0.14446444809436798
32枚のうち 16 枚が表の確率 = 0.13994993409141898
34枚のうち 17 枚が表の確率 = 0.13583375955931842
36枚のうち 18 枚が表の確率 = 0.13206059957155958
38枚のうち 19 枚が表の確率 = 0.1285853206354659
40枚のうち 20 枚が表の確率 = 0.12537068761957926
42枚のうち 21 枚が表の確率 = 0.12238567124768451
44枚のうち 22 枚が表の確率 = 0.11960417871932805
46枚のうち 23 枚が表の確率 = 0.11700408787760352
48枚のうち 24 枚が表の確率 = 0.11456650271348678
50枚のうち 25 枚が表の確率 = 0.11227517265921705

 つまり枚数が増えていくたびに、その半分が表になる確率はどんどん減っていく。

 余談だが、海外ではこういった「思い込みで 1/2 に見える確率」のようなイカサマ博打でカモられることもあるそうだ(日本でも江戸時代まではそういうことがあったらしい)。まぁ、賭け事なんて基本的には胴元が儲かる仕組み(だいたい胴元6:客4 以上の確率にしておけば、大数の法則で胴元は常に儲かる)になっているものなんだけどね。思い込みでダマされないようにしよう(笑)。


 最後に任意精度整数[BigInteger]とプリミティブ型での組合せの計算(long を double に書き換えたものも含む)での計算結果を比較してみよう。どのコードを使うかは求める答えや実行速度などに合わせて、自由に選べば良いと思う。

●long, double, BigInteger で組合せ計算した結果の比較
70C1
long = 70
double = 70
BigInteger = 70
70C2
long = 2415
double = 2415
BigInteger = 2415
70C3
long = 54740
double = 54740
BigInteger = 54740
70C4
long = 916895
double = 916895
BigInteger = 916895
70C5
long = 12103014
double = 12103014
BigInteger = 12103014
70C6
long = 131115985
double = 131115985
BigInteger = 131115985
70C7
long = 1198774720
double = 1198774720
BigInteger = 1198774720
70C8
long = 9440350920
double = 9440350920
BigInteger = 9440350920
70C9
long = 65033528560
double = 65033528560
BigInteger = 65033528560
70C10
long = 396704524216
double = 396704524216
BigInteger = 396704524216
70C11
long = 2163842859360
double = 2163842859360
BigInteger = 2163842859360
70C12
long = 10638894058520
double = 10638894058520
BigInteger = 10638894058520
70C13
long = 47465835030320
double = 47465835030320
BigInteger = 47465835030320
70C14
long = 193253756909160
double = 193253756909160
BigInteger = 193253756909160
70C15
long = 721480692460864
double = 721480692460864
BigInteger = 721480692460864
70C16
long = 2480089880334220
double = 2480089880334220
BigInteger = 2480089880334220
70C17
long = 7877932561061640
double = 7877932561061640
BigInteger = 7877932561061640
70C18
long = 23196134763125940
double = 23196134763125940
BigInteger = 23196134763125940
70C19
long = 63484158299081520
double = 63484158299081528 //←上から17桁以降は丸められる(0で埋まる[指数形式])
BigInteger = 63484158299081520
70C20
long = 161884603662657876
double = 161884603662657888
BigInteger = 161884603662657876
70C21
long = 385439532530137800
double = 385439532530137790
BigInteger = 385439532530137800
70C22
long = 19990591830327299 //←ここでオーバフロー
double = 858478958817125120
BigInteger = 858478958817125100
70C23
long = 41719495993726537
double = 1791608261879217410
BigInteger = 1791608261879217600
70C24
long = 81700679654381134
double = 3508566179513467400
BigInteger = 3508566179513467800
70C25
long = 150329250564061286
double = 6455761770304780300
BigInteger = 6455761770304780752
70C26
long = 260185241360875302
double = 11173433833219813000
BigInteger = 11173433833219812840
70C27
long = -259207164956705123
double = 18208558839321176000
BigInteger = 18208558839321176480
70C28
long = 260744142163258261
double = 27963143931814666000
BigInteger = 27963143931814663880
70C29
long = -258465175960438091
double = 40498346384007450000
BigInteger = 40498346384007444240
70C30
long = 261655728644386329
double = 55347740058143515000
BigInteger = 55347740058143507128
70C31
long = -257435965417228982
double = 71416438784701310000
BigInteger = 71416438784701299520
70C32
long = 262710669451175666
double = 87038784768854720000
BigInteger = 87038784768854708790
70C33
long = -256476928320147766
double = 100226479430802400000
BigInteger = 100226479430802391940
70C34
long = 263444050760708361
double = 109069992321755560000
BigInteger = 109069992321755544170
70C35
long = -256078807037830017
double = 112186277816662850000
BigInteger = 112186277816662845432


 組合せの計算に関しては、もう1つ「パスカルの三角形」を使った解法もあるので、ついでに覚えておくと色々役に立つかもしれない。


(関連記事)
【Java】パスカルの三角形を使って組合せを求める
【Java】順列の数を求める
【Java】プリミティブ型での階乗計算
【Java】1000000!のような巨大な階乗計算をする [任意精度整数(多倍長整数)]
【Java】任意精度整数(多倍長整数)で演算する
【Ruby】階乗・順列・組合せの個数を求める


category: Java

thread: プログラミング

janre: コンピュータ

tag: 算術関数  練習問題 
tb: 0   cm: --


プロフィール

Social

検索フォーム

全記事一覧

カテゴリ

ユーザータグ

最新記事

リンク

PR

▲ Pagetop