fc2ブログ
ヽ|∵|ゝ(Fantom) の 開発blog? ホーム »Java
カテゴリー「Java」の記事一覧

【Java】ジェネリクスで動的に配列を生成する  


 C# を使ってるとジェネリック型でも「T[] arr = new T[5];」のように簡単に配列を生成できるので、つい Java でもやろうとしてしまうね(笑)。残念ながら Java ではそういう書き方はできないのと、プリミティブ型はジェネリクスには当てはめられない欠点があるので、専用メソッドを作って置くと便利だ。型変換については必要なだけオーバーロードでも作るしかないが、とりあえず簡単な例を挙げておこう。応用すれば色々役に立つだろう。


(※) Java openjdk 9 (build 9+181)[paiza.io] / Android Studio 3.1.2 [Java 1.7] / Windows10(x64) で確認



■ジェネリクスで動的に配列を生成する

●ジェネリクスで配列を動的に生成する(関数定義)
import java.lang.reflect.Array;

//1次元配列を生成する
@SuppressWarnings("unchecked")
public static final <T> T[] newArray(final int length, final T[] type) {
return (T[]) Array.newInstance(type.getClass().getComponentType(), length);
}

//2次元配列を生成する
@SuppressWarnings("unchecked")
public static final <T> T[][] newArray(final int y, final int x, final T[] type) {
return (T[][]) Array.newInstance(type.getClass().getComponentType(), y, x);
}

//3次元配列を生成する
@SuppressWarnings("unchecked")
public static final <T> T[][][] newArray(final int z, final int y, final int x, final T[] type) {
return (T[][][]) Array.newInstance(type.getClass().getComponentType(), z, y, x);
}

●使用例(メインでのコードなど)
String[] arr1 = newArray(5, new String[0]);  //new String[5] と同じ
System.out.println("length = " + arr1.length);
System.out.println();

String[][] arr2 = newArray(2, 5, new String[0]); //new String[2][5] と同じ
System.out.println("y length = " + arr2.length);
System.out.println("x length = " + arr2[0].length);
System.out.println();

String[][][] arr3 = newArray(3, 2, 5, new String[0]); //new String[3][2][5] と同じ
System.out.println("z length = " + arr3.length);
System.out.println("y length = " + arr3[0].length);
System.out.println("x length = " + arr3[0][0].length);
System.out.println();

length = 5

y length = 2
x length = 5

z length = 3
y length = 2
x length = 5

 引数の数値は次元数だが、type の引数に「new String[0]」を与えているのに注意しよう。これは配列そのものを渡しているわけでなく、配列の要素の型を取得して、ジェネリクスに当てはめるために使用するものだ。この方法はジェネリクスの型を取得する裏技的な手法として Java では良く使われる(単一で型を取得するメソッドはないが、配列で要素の型を取得するメソッドはあるため(Class.getClass().getComponentType())。「T...」みたいな可変長引数でも良い)。List を配列に変換するメソッド「List.toArray(new String[0])」を使ったことがあるかも知れないが、それと用法は同じだ。なので要素の数は0で良い。



■Integer[] → int[] 変換

 次に int[] を生成する方法を考えてみよう。残念ながら Java では元々プリミティブ型(int, float, boolean など)はジェネリクスに使えないという欠点がある。なのでとりあえずラッパークラスの Integer の配列で生成し、int[] 型に変換するという方法をとれば良い。あまりスマートではない気もするが、簡単なのでロジックで悩むことは無いだろう。

●Integer[] → int[] に変換する(関数定義)
//Integer[] → int[] に変換する
public static final int[] toIntArray(final Integer[] src) {
final int len = src.length;
final int[] dst = new int[len];
for (int i = 0; i < len; i++)
dst[i] = (src[i] == null) ? 0 : src[i];
return dst;
}

●使用例(メインでのコードなど)
//Integer[] で生成(※要素が null なので使う際は注意)
Integer[] tmp = newArray(5, new Integer[0]);
for (int i = 0; i < tmp.length; i++)
tmp[i] = i; //連番で要素を初期化

//Integer[] → int[] に変換
int[] arr = toIntArray(tmp);
System.out.println("length = " + arr.length);

//要素を結合して表示
String s = "";
for (int e : arr) {
if (s.length() > 0)
s += ", ";
s += e;
}
System.out.println(s);

length = 5
0, 1, 2, 3, 4

 またはラムダ式の使える環境(Java 1.8 / Android 7.0[API 24])なら「Integer[] → int[] 変換」を以下のようにも書ける。

●ラムダ式(steam)を使った「Integer[] → int[] 変換」
import java.util.Arrays;

//Integer[] → int[] に変換
int[] arr = Arrays.stream(tmp).mapToInt(e -> e).toArray(); //Java 1.8 / Android 7.0[API 24]

 ラムダ式は Java 1.8 以降、Arrays.steam() は Android 7.0[API 24] 以降なので、スマホの古い機種も対応させたいなら使わない方が良いかも知れない(Android って現時点ではまだ Java 1.7 が多いので)。Java のジェネリクスはプリミティブ型が使えないのが痛いが(だから配列系のメソッドにはオーバーロードが鬼のようにある(笑))、他のプリミティブ型も同じようにラッパークラスを経由すれば、動的生成も可能だ。簡単なのでやってみると良いだろう。


(関連記事)
【Java】配列, リスト(List), 連想配列(Map) の初期化
【Java】配列要素の反転(reverse)
【Java】2次元配列のソート


スポンサーサイト



category: Java

thread: プログラミング

janre: コンピュータ

tag: Javaライブラリ 
tb: 0   cm: --

【Java】数値の桁数を調べる(べき乗の桁数・べき乗のべき乗の桁を調べる)  


 数値の桁数を知りたいシーンはわりと良くあるもので、調べてみると「10で割っていく方法」や「文字列に変換する方法」「対数を使う方法」「配列にあらかじめ各桁の最大(9999のような)を入れておいて比較する方法」などあるようだが、実行速度は言語によってまちまちのようなので、簡単なアルゴリズムを色々知っておいても損はないだろうと試してみた。




■10で割っていく方法

 まずは一番オーソドックスな方法(?)でやってみよう。これは以前書いた「基数変換のアルゴリズム」と同じと言っても良い。基数変換では具体的な桁の値をリストで返しているが、それをカウントに変えるだけだ。故に基数変換の戻値のリストのサイズ(List.size())でも同じとなる。まずは基数変換のコードの一部を抜き出して少しばかり改造してみよう。

●10で割って桁数をカウントする
/**<h1>整数の桁数を求める</h1>
* <p>0 となるまで、10 で割り続ける方法。</p>
* @param value : 調べる値
* @return<b>int</b> : 桁数 (負の値は '-' を除いた桁数)
*/
public static final int length(long value) {
int cnt = 0;
do {
cnt++;
value /= 10;
} while (value != 0);
return cnt;
}


//メインでは...
int n = 1234; //4
//int n = -4321; //4
//int n = 0; //1
System.out.println(length(n));

4

 負の値のとき、マイナス記号 '-' を除外してカウントしているが、負の割り算は言語によって仕様が違うので、他の言語に移植するときは注意しよう(Ruby などは上手くいかない)。心配なら関数のはじめに絶対値(Math.abs())などで正の値にしておくと良い(Java では無くても良い)。調べる値が 0 のとき、1 桁としたいので、すぐにインクリメントしているが、while に書きなおしたいときは、はじめからカウントを 1 にしておけば良いだろう。上記の例は基数変換のコードをもとにしているが、以下のように while でシンプルに書くこともできる。

●上記の例を while に書きなおしたもの
/**<h1>整数の桁数を求める</h1>
* <p>0 となるまで、10 で割り続ける方法。</p>
* @param value : 調べる値
* @return<b>int</b> : 桁数 (負の値は '-' を除いた桁数)
*/
public static final int length(long value) {
int cnt = 1; //0 にも対応
while ((value /= 10) != 0) {
cnt++;
}
return cnt;
}


//メインでは...
int n = 1234; //4
//int n = -4321; //4
//int n = 0; //1
System.out.println(length(n));

4

 内容は全く同じで、実行速度も変わらないので、好みの書き方で選べば良い。




■文字列に変換する方法

 「10で割る方法」と同じようによく使われるのは、「文字列に変換する方法」だ。ただし、Java の場合は文字列変換は負荷が高いので、実行速度は遅くなる。しかし、Ruby をはじめ、PHP や Perl のようなスクリプト言語は文字列変換のスピードは非常に速いので、むしろ「10で割る方法」の方が遅いらしい。参考までに色々な言語でやっているページを載せておくが、私が実際に Ruby で試したところ、「loop do~break if~」の代わりに「while」に書きなおしたら、ここまでの差は出なかった(それでも文字列変換の方がわずかに速い)。他の言語から移植するときのためにも、この手のやり方には慣れておいた方が良いだろう。

(参考)
整数の桁数を求める方法

●文字列に変換し、長さを取得する
//int n = 1234;     //4
int n = -4321; //4
//int n = 0; //1

System.out.println(String.valueOf(Math.abs(n)).length());

4

 なお、値が正の数だけとわかっていれば、Math.abs() は必要ない。




■べき乗(xk)の桁数を調べる

 次に、普通に計算したらオーバーフローするような「べき乗の値の桁を調べる方法」をやってみよう。それは(常用)対数関数Math.log10())を使う方法だ。ちなみにこの方法で int や long の桁数を調べることもできるが、log() は負荷が高い関数なので、「10で割る方法」の方がよほど速い。なので、整数型では求められない値など、使い分けが必要になるだろう。

 まずは簡単に対数関数の使い方をおさらいしてみよう。といっても高校の数学の教科書そのままの内容だ(笑)。また対数関数指数関数は裏返しの関数なので、ついでにおさらいしておくとなお良い。以下の問題の解き方を丸暗記してしまえば、わりと簡単に理解できるハズだ。

【問】230は何桁の数字か。log102 = 0.3010 を用いて調べよ。

【解】
log102 = 0.3010 だから、2 = 100.3010  [指数への変換:logap = q ⇔ p = aq]
230 = (100.3010)30
  = 100.3010×30  [指数法則:(ap)q = ap×q]
  = 109.03
よって、109 < 230 < 1010
したがって、230 は 10桁の数字である。  [109 = 1,000,000,000 (頭 1 と 0 が 9 個で 10桁)]


 さて、上の解法を理解したら、いよいよ本題のべき乗(xk)の桁数を求めてみよう。例で使ってるのは底が10である常用対数の使い方であるが、それは言い換えればある数値 n を 10q に置き換える関数であり、桁を表していると言える。つまり、ここでは具体的な値ではなく、桁数を求めたいわけだから、10 の指数部分 q を求めれば良いわけである。

 上の例を用いれば桁数は 230 = 109.03 の 9.03 部分のことで、これは log102×30 であり、この値の小数部分を切り捨て、1 を足したものである(頭 1 と 0 が 9 個で 10桁)。この 2 を x、30乗を k と置けば、以下のようになっているとわかるだろう。

230 = 100.3010×30 = 10log102×30

(2 を x、30乗を k と置くと)
xk = 10log10x × k  … ① (←以降、この式を公式のように扱う)

(桁数は指数部分のことだから)
桁数 = floor(log10x × k) + 1

●べき乗(xk)の桁数を log10() により求める
/**<h1>x<sup>n</sup> の桁数を求める</h1>
* <p>log10() により求める。底が負のときは NaN となるので注意。</p>
* @param x : 底 (> 0)
* @param k : べき指数 (>= 0)
* @return<b>double</b> : 桁数
*/
public static final double powLength(final double x, final double k) {
return Math.floor(Math.log10(x) * k) + 1;
}


//メインでは...
int x = 2, k = 30; //10.0
//int x = 123456789, k = 99999; //809144.0
System.out.println(powLength(x, k));

10.0

 log() の引数は負の値だと NaN が返される注意しよう。必要ならエラーを投げるのも良い。引数の値によっては非常に大きな桁数になるので double のまま返しているが、値の範囲がわかっていれば整数型にキャストするのも良いだろう。Java の場合、long は 19桁の値、double は指数表現を含めて 309桁の値まで表せるが、べき乗はあっという間にオーバーフロー(Infinity)することがあるので、使う際にはあらかじめ確認しておいた方が良いだろう。




■べき乗のべき乗(xyz)の桁数を求める

 次に更にべき乗をかけた「べき乗のべき乗」の桁を調べてみよう。考え方は「べき乗(xk)の桁数を調べる」をそのままで、少しばかり拡張したものになる。その式を考えてみよう。

(xyz の指数部分 yz を log10 を使って変換すると [①を参照])
yz = 10log10y × z

(xyz の指数部分 yz を k と置き、同じように変換すると [①を参照])
xyz = xk = 10log10x × k

(k を元に戻すと)
xyz = xk = 10log10x × 10log10y × z

(桁数は指数部分のことだから)
桁数 = floor(log10x × 10log10y × z) + 1

●べき乗のべき乗(xyz)の桁数を log10() により求める
/**<h1>x<sup>y<sup>z</sup></sup> (べき乗のべき乗)の桁数を求める</h1>
* <p>log10() により求める。底が負のときは NaN となるので注意。</p>
* @param x : 底 (> 0)
* @param y : べき指数1 (>= 0)
* @param z : べき指数2 (>= 0)
* @return<b>double</b> : 桁数 (オーバーフロー: Infinity に注意)
*/
public static final double powpowLength(final double x, final double y, final double z) {
final double p = Math.log10(y) * z; //y^z
final double q = Math.log10(x); //x^p
return Math.floor(q * Math.pow(10, p)) + 1;
}


//メインでは...
int x = 2, y = 3, z = 5; //74.0
//int x = 5, y = 6, z = 7; //195667.0
System.out.println(powpowLength(x, y, z));

74.0




■2つのべき乗のべき乗(xyz, abc)の大きさを比較する

 ここでは2つのべき乗のべき乗(xyz, abc)の大小関係を比較することを考えてみよう。もちろん前述した「べき乗のべき乗(xyz)の桁数を求める」で比較しても構わないが、なぜわざわざ比較だけをピックアップしたかというと、引数によっては「べき乗のべき乗」では桁数が非常に大きくなりオーバーフロー(Infinity)しかねないので、対数のまま利用した方が良いからだ。もう一度「べき乗のべき乗(xyz)の桁数を求める」のコードを見てみよう。桁数が大きくなる原因に「Math.pow(10, p)」という部分がある。ここを変形してオーバーフローを防ぐ方法を考えてみよう。

(桁数を表す指数部分を抜き出すと)
yz = log10x × 10log10y × z

log10x = m と置き、今までと同じように[①を参照] 10q に変形すると、m = 10log10m [m1 = 10log10m × 1

(m を元に戻すと)
10log10m = 10log10(log10x) (= log10x)

(yzの式に代入して)
yz = 10log10(log10x) × 10log10y × z = 10log10(log10x) + log10y × z [指数法則:ap×aq = ap+q]

(xyzの式に戻すと)
xyz = 1010log10(log10x) + zlog10y

(つまり、一番上の指数だけを比較すれば大小関係がわかる)
log10(log10x) + zlog10y

●べき乗のべき乗(xyz)の比較をする
/**<h1>x<sup>y<sup>z</sup></sup> (べき乗のべき乗) の log10() を返す</h1>
* <p>この値で比較すれば、大小関係がわかる。底が負のときは NaN となるので注意。</p>
* @param x : 底 (> 0)
* @param y : べき指数1 (>= 0)
* @param z : べき指数2 (>= 0)
* @return<b>double</b> : 底を10<sup>10</sup>に合わせた log10() を返す
*/
public static final double powpowLog10(final double x, final double y, final double z) {
final double p = Math.log10(y) * z; //y^z
final double q = Math.log10(Math.log10(x)); //10^q = log10 x
return p + q;
}

/**<h1>x<sup>y<sup>z</sup></sup> (べき乗のべき乗) 同士の大小の比較</h1>
* <p>比較の関係値を返す(-1, 0, 1 のみ)。差分ではない。</p>
* @param x1 : 底 a (> 0)
* @param y1 : べき指数1 a (>= 0)
* @param z1 : べき指数2 a (>= 0)
* @param x2 : 底 b (> 0)
* @param y2 : べき指数1 b (>= 0)
* @param z2 : べき指数2 b (>= 0)
* @return<b>int</b> : -1 : a < b / 1 : a > b / 0 : a == b
*/
public static final int powpowCompare(final double x1, final double y1, final double z1,
final double x2, final double y2, final double z2) {

final double a = powpowLog10(x1, y1, z1);
final double b = powpowLog10(x2, y2, z2);
if (a < b) {
return -1;
} else if (a > b) {
return 1;
} else {
return 0;
}
}


//メインでは...
int x1 = 123456, y1 = 9876543, z1 = 24356891; //Infinity
int x2 = 333333, y2 = 5555555, z2 = 22222222; //Infinity
System.out.println(powpowLength(x1, y1, z1));
System.out.println(powpowLog10(x1, y1, z1));
System.out.println(powpowLength(x2, y2, z2));
System.out.println(powpowLog10(x2, y2, z2));
System.out.println(powpowCompare(x1, y1, z1, x2, y2, z2));

Infinity
1.7036683127845693E8
Infinity
1.4988283149816477E8
1

 この方法で行けば、309桁(Double.MAX_VALUE = 1.7976931348623157E308) までは比較できる

 また、比較だけなら、Math.log10() の代わりに Math.log() を使っても良いが、値が大きくなる可能性があるので、Math.log10() を使った方が無難かも知れない(実行速度も変わらない)。もちろん引数が負の値のときは NaN となるので注意しよう。戻値は比較の関係値を返す -1, 0, 1 のみにしているが、必要なら (a - b) の差分を返すように改造しても良いだろう(ただし、型に注意。比較関数は大抵 int になっているので)。

 ちなみに同じように「べき乗のべき乗のべき乗」も考えられるが、log() を3回以上入れ子にすると負の値となってしまうので、上手く行かないようだ。式としては以下のようになる。

wxyz = 101010log10(log10(log10w)) + log10(log10x) + zlog10y

(一番上の指数だけを抜き出せば)
log10(log10(log10w)) + log10(log10x) + zlog10y

 式は無理矢理なので、正しい書き方とかは勘弁して欲しい(笑)。級数(数列の和)っぽく書けばもっとカッコイイのかも知れないが、この書き方なら高校数学の範疇なので理解できるハズだ(無理矢理感は別として)。でもまぁ、log() でもこれ以上は無理そうなので、他の方法を研究してみるのも面白いだろう。




■任意精度を使う

 最後に任意精度実数を使って桁を調べてみよう。考え方としては「文字列に変換する方法」と同じと思って良い。具体的な値を返すので、桁が大きくなるほど重くなるが、検算したいときには役に立つこともあるだろう(といっても100万桁ぐらいが限界でそれ以上はフリーズしたようになる)。なお、BigInteger を用いても構わないが、桁を調べるだけなら、文字列に変換する(toString())する必要はないので、整数部分の桁を返す BigDecimal.precision() を使った方が速い

●任意精度実数で具体的な値と桁数を調べる
int x = 2, y = 3, z = 5;
int a = (int)Math.pow(y, z);
BigDecimal b = new BigDecimal(x).pow(a);
System.out.println(b.toString());
System.out.println(b.precision());
System.out.println(powpowLength(x, y, z));

14134776518227074636666380005943348126619871175004951664972849610340958208
74
74.0

 上の例はべき乗のべき乗(235)と同じものだ。ただし、BigDecimal.pow() の引数は int 型にしか対応していない(と言っても、long でできたとして重すぎてフリーズしたようになると思うが)。具体的な値を知りたければ文字列化(toString())すれば良いが、桁が大きいときは負荷が高いので、あらかじめ値の範囲を調べた上で使った方が良いだろう。オンラインジャッジでの問題などでは、任意精度ではほとんどの場合タイムアウトすると思うので、桁数だけなら、前述の「powLength()」(べき乗の桁数)や「powpowLength()」(べき乗のべき乗の桁数)を使った方が速い。





(関連記事)
【Java】n進数変換(基数変換)・桁の分割をする
【Java】べき乗計算を高速化する(繰り返し二乗法)
【Java】任意精度整数(多倍長整数)で演算する


category: Java

thread: プログラミング

janre: コンピュータ

tag: 算術関数  べき乗  アルゴリズム 
tb: 0   cm: --

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


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

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



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

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

(漸化式)
F0 = 0,
F1 = 1,
Fn = Fn-1 + Fn-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 でオーバーフローになるので気をつけよう。


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

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


●ビネの公式でフィボナッチ数を求める
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() でも整数型のキャストでも良いだろう(正の整数の場合)。よくわからない記号は「数学記号の表」などで調べると良い。

 ビネの公式での誤差の出る範囲を抑える1つの方法は、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 で計算した結果をダンプしてみよう。オーバーフローや誤差を確認しておけば、どれを使うかの目安になるかも知れない。値を確認するには以下のようなサイトを参考にすれば良い。

(参考)[数列]
100番目までのフィボナッチ数列
1番目から1000番目までのフィボナッチ数列

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】3つ以上の最大公約数を求める(ユークリッドの互除法②)  


 「ユークリッドの互除法」や「拡張ユークリッドの互除法」は「自然数 a, b 2つの最大公約数を求める」ものであるが、では3つ以上の複数の最大公約数(a, b, c, d,...)はどうやって求めるのだろうか。実はこれも英語版の Wikipedia に書いてある。説明を翻訳しなくても「gcd(a, b, c) = gcd(gcd(a, b), c)」の式を見れば十分わかるだろう。簡単に言えば「a, b の最大公約数を求め、更にその答えと c の最大公約数を求める」ということになる。更に複数のパラメタがあるときは、それを繰り返す感じになる。つまりは入れ子にすれば良いだけである。


 といっても、3つ4つ…と増やした場合は表記が面倒になるので、いっそのこと専用の関数を作ってしまおう。Java には可変引数(型...)があるので、それを使えば簡単に書ける。

●3つ以上の最大公約数を求める[ユークリッドの互除法②]
/**<h1>3つ以上の最大公約数を求める[gcd(a, b, c,...)]</h1>
* <p>ユークリッドの互除法の繰り返しで求める。</p>
* @param param : 数値1, 2, 3,... (a, b, c,...) (>0) [引数は3個以上]
* @return<b>int</b> : 最大公約数(なし=1[互いに素])
*/
public static final int gcd(final int... param) {
final int len = param.length;
int g = gcd(param[0], param[1]); //gcd(a, b) は以前作ったもの
for (int i = 1; i < len - 1; i++) {
g = gcd(g, param[i + 1]); //gcd(a, b) は以前作ったもの
}
return g;
}


//メインでは...
int a = 12;
int b = 36;
int c = 60;
int d = 96;

int g1 = gcd(a, b);
int g2 = gcd(g1, c);
int g3 = gcd(g2, d);
System.out.println("g3 = " + g3);

int g = gcd(a, b, c, d);
System.out.println("g = " + g);
System.out.println(g == g3);

g3 = 12
g = 12
true

 gcd(a, b) は以前作ったものをそのまま使って欲しい。入れ子はループで書いてあるが、いかにも再帰で書けそうな感じなので、試してみても良いかも知れない。ここでは後述する「3つ以上のパラメタに対応させた拡張ユークリッド互除法」と比較するためにもループで書いてみる。

 また、引数が1つのときは不正となるので気をつけよう。必要があれば例外処理を入れても良いかも知れない(引数が1つ=param.length(len)が1のときはエラーとするなど)。gcd() を同名でオーバーロードした場合は引数2つの場合は gcd(a, b) になり、3つ以上の場合は gcd(...) の方にディスパッチされる。ちなみに可変引数はメソッド内では配列のように扱われる。紛らわしいと感じるなら別名で関数定義しても良いだろう。

 それではもう1つ、同じように「拡張ユークリッドの互除法」を入れ子にして、複数のパラメタに対応させてみよう。

●3つ以上のパラメタに対応させた拡張ユークリッドの互除法
/**<h1>3つ以上のパラメタに対応させた拡張ユークリッド互除法[ax + by + cz +... = gcd(a, b, c,...)]</h1>
* <p>拡張ユークリッドの互除法の繰り返しで求める。</p>
* @param param : 数値 1, 2, 3,... (a, b, c,...) (>0) [引数は3個以上]
* @return<b>int[]</b> : {[0]:gcd, [1]:x, [2]:y, [3]z,...} 最大公約数(なし=1 [互いに素])と解 x, y, z...
*/
public static final int[] extgcd(final int... param) {
final int len = param.length;
final int[][] g = new int[len - 1][];
g[0] = extgcd(param[0], param[1]); //extgcd(a, b) は以前作ったもの
for (int i = 1; i < len - 1; i++) {
g[i] = extgcd(g[i - 1][0], param[i + 1]); //extgcd(a, b) は以前作ったもの
}

//解の係数を展開する
final int[] res = new int[len + 1];
res[len] = g[len - 2][2];
res[0] = g[len - 2][0]; //gcd

int k = g[len - 2][1];
for (int i = len - 3; i >= 0; i--) {
res[i + 2] = g[i][2] * k;
k *= g[i][1];
}
res[1] = k;
return res;
}


//メインでは...
int a = 12;
int b = 16;
int c = 18;
int d = 21;

int[] g1 = extgcd(a, b);
System.out.println("g1 = " + g1[0] + ", x = " + g1[1] + ", y = " + g1[2]);
int[] g2 = extgcd(g1[0], c);
System.out.println("g2 = " + g2[0] + ", x = " + g2[1] + ", y = " + g2[2]);
int[] g3 = extgcd(g2[0], d);
System.out.println("g3 = " + g3[0] + ", x = " + g3[1] + ", y = " + g3[2]);

int[] g = extgcd(a, b, c, d);
System.out.println("g = " + g[0] + ", w = " + g[1] + ", x = " + g[2] + ", y = " + g[3] + ", z = " + g[4]);

g1 = 4, x = -1, y = 1
g2 = 2, x = -4, y = 1
g3 = 1, x = -10, y = 1
g = 1, w = -40, x = 40, y = -10, z = 1

 extgcd(a, b) は以前作ったものをそのまま使って欲しい。これも gcd(...) と基本的には変わらないので、引数が1つのときはエラーとなるので注意しよう。例外処理を追加しても良い。

 拡張ユークリッドの互除法を繰り返している部分は gcd(...) はほぼ同じと言っても良いが、解も入れ子になっているので、それらを展開するために配列に1回ごとの extgcd(a, b) の答えを記録し、最後に計算している。コードが少し難しく感じるかも知れないが、元々は以下のように、1つずつの式だと思えばそれほどでないことがわかるだろう。

a = 12, b = 16, c = 18, d = 21 のとき、
g1 = extgcd(a, b) → g1 = 4, x = -1, y = 1 → 4 = 12 * -1 + 16 * 1
g2 = extgcd(g1, c) → g2 = 2, x = -4, y = 1 → 2 = 4 * -4 + 18 * 1
g3 = extgcd(g2, d) → g3 = 1, x = -10, y = 1 → 1 = 2 * -10 + 21 * 1

g = extgcd(extgcd(extgcd(a, b), c), d) のように入れ子になっているので、
1 = ((12 * -1 + 16 * 1) * -4 + 18 * 1) * -10 + 21 * 1 (式を代入する)
 = 12 * -40 + 16 * 40 + 18 * -10 + 21 * 1 (展開する)
g = a * w + b * x + c * y + d * z (対応させる)

 使う際に少しだけ気をつけないといけないのは、係数(a, b, c, d)の値の順番を入れ替えると解(w, x, y, z)は変わる点だ。というのは a, b で求めた解 w, x を用いて更に y を求め、更に z を求め…という感じで計算しているので、1つ前の答えによって後の答えが変わるからだ。もちろんどの式も検算すれば正しいことがわかる。「いちばん小さい値(全ての解の絶対値の和とか)を求めよ」などのときには、最小となるパターンを別途見つけなくてはならない。例えば、上の例の係数を順列で並び変えてみると、

a = 12, b = 16, c = 18, d = 21, g = 1, w = -40, x = 40, y = -10, z = 1
a = 12, b = 16, c = 21, d = 18, g = 1, w = 5, x = -5, y = 1, z = 0
a = 12, b = 18, c = 16, d = 21, g = 1, w = 30, x = -30, y = 10, z = 1
a = 12, b = 18, c = 21, d = 16, g = 1, w = -15, x = 15, y = -5, z = 1
a = 12, b = 21, c = 16, d = 18, g = 1, w = -10, x = 5, y = 1, z = 0
a = 12, b = 21, c = 18, d = 16, g = 1, w = -10, x = 5, y = 0, z = 1
a = 16, b = 12, c = 18, d = 21, g = 1, w = 40, x = -40, y = -10, z = 1
a = 16, b = 12, c = 21, d = 18, g = 1, w = -5, x = 5, y = 1, z = 0
a = 16, b = 18, c = 12, d = 21, g = 1, w = 10, x = -10, y = 0, z = 1
a = 16, b = 18, c = 21, d = 12, g = 1, w = 10, x = -10, y = 1, z = 0
a = 16, b = 21, c = 12, d = 18, g = 1, w = 4, x = -3, y = 0, z = 0
a = 16, b = 21, c = 18, d = 12, g = 1, w = 4, x = -3, y = 0, z = 0
a = 18, b = 12, c = 16, d = 21, g = 1, w = -30, x = 30, y = 10, z = 1
a = 18, b = 12, c = 21, d = 16, g = 1, w = 15, x = -15, y = -5, z = 1
a = 18, b = 16, c = 12, d = 21, g = 1, w = -10, x = 10, y = 0, z = 1
a = 18, b = 16, c = 21, d = 12, g = 1, w = -10, x = 10, y = 1, z = 0
a = 18, b = 21, c = 12, d = 16, g = 1, w = 5, x = -5, y = 0, z = 1
a = 18, b = 21, c = 16, d = 12, g = 1, w = 5, x = -5, y = 1, z = 0
a = 21, b = 12, c = 16, d = 18, g = 1, w = 5, x = -10, y = 1, z = 0
a = 21, b = 12, c = 18, d = 16, g = 1, w = 5, x = -10, y = 0, z = 1
a = 21, b = 16, c = 12, d = 18, g = 1, w = -3, x = 4, y = 0, z = 0
a = 21, b = 16, c = 18, d = 12, g = 1, w = -3, x = 4, y = 0, z = 0
a = 21, b = 18, c = 12, d = 16, g = 1, w = -5, x = 5, y = 0, z = 1
a = 21, b = 18, c = 16, d = 12, g = 1, w = -5, x = 5, y = 1, z = 0

min = 7
a = 16, b = 21, c = 12, d = 18, g = 1, w = 4, x = -3, y = 0, z = 0
a = 16, b = 21, c = 18, d = 12, g = 1, w = 4, x = -3, y = 0, z = 0
a = 21, b = 16, c = 12, d = 18, g = 1, w = -3, x = 4, y = 0, z = 0
a = 21, b = 16, c = 18, d = 12, g = 1, w = -3, x = 4, y = 0, z = 0

のようになり、最小となる解の絶対値の総和は7となる(解に0を許す場合)。

 再帰でも書いてもいいかも知れないが、特に拡張ユークリッド互除法の方はコードが複雑化し兼ねないので気をつけよう。再帰よりループの方があとから改造しやすい利点もあるしね。





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


category: Java

thread: プログラミング

janre: コンピュータ

tag: 算術関数  約数 
tb: 0   cm: --

【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: --


プロフィール

Social

検索フォーム

全記事一覧

カテゴリ

ユーザータグ

最新記事

リンク

PR

▲ Pagetop