【Java】３つ以上の最大公約数を求める（ユークリッドの互除法②）

2016/06/20 Mon [edit]

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

■３つ以上のユークリッドの互除法
■３つ以上の拡張ユークリッドの互除法

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

●３つ以上の最大公約数を求める[ユークリッドの互除法②]

/**<h1>３つ以上の最大公約数を求める[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) は以前作ったものをそのまま使って欲しい。入れ子はループで書いてあるが、いかにも再帰で書けそうな感じなので、試してみても良いかも知れない。ここでは後述する「３つ以上のパラメタに対応させた拡張ユークリッド互除法」と比較するためにもループで書いてみる。

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

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

●３つ以上のパラメタに対応させた拡張ユークリッドの互除法

/**<h1>３つ以上のパラメタに対応させた拡張ユークリッド互除法[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(...) と基本的には変わらないので、引数が１つのときはエラーとなるので注意しよう。例外処理を追加しても良い。

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

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 を求め…という感じで計算しているので、１つ前の答えによって後の答えが変わるからだ。もちろんどの式も検算すれば正しいことがわかる。「いちばん小さい値（全ての解の絶対値の和とか）を求めよ」などのときには、最小となるパターンを別途見つけなくてはならない。例えば、上の例の係数を順列で並び変えてみると、

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

のようになり、最小となる解の絶対値の総和は７となる（解に０を許す場合）。

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