- できるだけ毎週土曜日更新．次回更新は2月16日．
- ネタ切れ深刻なので記事リクエスト募集中．twitterないしはコメント欄まで．
- 貼っているコードはあんまり検証していないので，信頼性はアレです．アレ．
- twitter: @tmaehara フォローはお気軽に．
2013-01-26
Euclidean algorithm (Blankinship algorithm)
a, b を互いに素な整数としたとき，方程式 a x + b y = 1 の整数解を求める手続きは拡張ユークリッド互除法（extended Euclidean algorithm）として知られています．今回はこれに関連する話題を紹介します．
https://topcoder-g-hatena-ne-jp.jag-icpc.org/iwiwi/20130105/1357363348 で Mi_Sawa さんが紹介している非再帰版のユークリッド互除法（以下の手続き）を導出するところからはじめます．
LL extGCD(LL a, LL b, LL &x, LL &y) {
for (LL u = y = 1, v = x = 0; a; ) {
LL q = b / a;
swap(x -= q * u, u);
swap(y -= q * v, v);
swap(b -= q * a, a);
}
return b;
}
まず次の行列を考えます．
A' = |a 1 0|
|b 0 1|
（プライムは後の都合による）．すると，x = [-1, a, b] に対して A' x = 0 が成立します．そこで，この等式を崩さないように A' を行基本変形していきます．整数行列に対する行基本変形は
- ある行に別の行の整数倍を加える
- ある行の符号をかえる
- ある行と別の行を入れ替える
という操作の繰り返しでした．これらは全て行列の左から 2×2 行列 S をかけることで実現できます．A' x = 0 の左から S をかけると S A' x = 0 なので，行基本変形は「x を解とする」ことを保ちます．また，基本変形がすべて可逆なので S は正則であり，S A' x = 0 から A' x = 0 が導けることも注意しておきます．
さて，行基本変形の目的地として，一番左の列が簡単な形になるようにしましょう．基本変形の操作でユークリッドの互除法がシミュレートできることに注意すると，基本変形で
S A' = |g x y|
|0 u v|
までいけることがわかります．ここで g = gcd(a,b) です．これが [-1,a,b] を解にもつ，ということを書くと
a x + b y = g a u + b v = 0
となって不定方程式の解がもとまっていることがわかります．以上の手続きをコードにすると冒頭のものになり，非再帰の拡張互除法が導出できました．
...
Blankinship's algorithm
上の形でのユークリッドの互除法の理解は [Blankinship 1963] によります．一般に，不定方程式
a_1 x_1 + ... + a_n x_n = gcd(a_1, ..., a_n)
を解くためには行列 A' = [A, I] （A は a が縦に並んだ行列，右は辻褄の合う単位行列）について基本変形を行えばよい，ということを Blankinship は指摘しました．この意味で冒頭の非再帰版の拡張互除法は Blankinship's algorithm とも呼ばれます．
なお，この手続きによって実は不定方程式の一般解も求まっている，ということが [Morito-Salkin 1979] によって指摘されています．上の２変数の例でいえば，すべての解は [x,y] + k [u,v] で書ける，ということですね．
typedef long long LL; LL extGCD(LL a[], LL x[], int n) { LL A[n][n+1]; memset(A, 0, sizeof(A)); for (int i = 0; i < n; ++i) { A[i][n] = a[i]; A[i][i] = 1; } while (1) { int k = -1; // k = nonzero argmin A[*][n] for (int i = 0; i < n; ++i) if (A[i][n]) if (k == -1 || abs(A[k][n]) > abs(A[i][n])) k = i; bool fin = true; for (int i = 0; i < n; ++i) { if (i == k || A[i][n] == 0) continue; fin = false; LL q = A[i][n] / A[k][n]; for (int j = 0; j <= n; ++j) A[i][j] -= q * A[k][j]; } if (fin) { for (int j = 0; j < n; ++j) // A[k] = x; a x = g (particular solution) x[j] = A[k][j]; // A[!k] = u; a u = 0 (homegeneous solution) return A[k][n]; } } } int main() { int n = 3; LL x[n], a[n] = {12,15,10}; LL g = extGCD(a, x, n); for (int i = 0; i < n; ++i) cout << x[i] << " "; // -2 1 1 }
線形方程式の求解として見る
※コードここまで
上の手続きの別の見方として，線型方程式 [x,y] A = c を解く問題だと考えます．右辺は拡張互除法では g でしたが，一般的に c にしておきます．普通の線型代数で方程式 x A = c を解こうと思ったら，A を適当に基本変形してくことになるので，それに倣って基本変形していきます．ただし，今 A はサイズ的に左からの積しか意味が無いので，左から操作します．左からの積は行基本変形に対応していたことに注意すると
A = |a|
|b|
に適当に左から行列 S をかける，ということなので
S A = |g|
|0|
となって，x S^{-1} = [s,t] とおけば結局 x A = s g = c となって，c が g で割り切れるときのみ整数解をもち，「s = 1/g, t = 任意」がその解全てであることがわかります．[s,t] から x を求めるには [x,y] = S [s,t] なので S が求まれば十分です．これは学部の線形代数の講義を思い出すと，A を A' = [A I] と拡張しておいて A 側に基本変形を行なって I に伝播させればよい，ということでした．この手続きは Blankinship の手続きと一致しています．
※学部の線型代数でやることは，行列 A の逆行列を求めるために [A I] を考えて基本変形で [I S] へと変形すると S = A^{-1} になる，というものだったはずです．今は A の逆行列まではいけないので「S A = 簡単な形」とする S をあらわに求める手続きになっています．
単因子標準形
さらに拡張します．A を整数係数の行列，b を整数ベクトルとして，線型方程式 x A = b を満たす x を求める問題を考えます．今度は A は行基本変形も列基本変形も意味があるので両方やります．どこまで A がいけるか，ということに答えるものが単因子標準形です．
定理． 整数を成分とする行列 A にたいし，あるユニモジュラ行列 S, T が存在して
- S A T = diag(g[1], ..., g[r]) （サイズのあわないところはゼロ）
- g[1] = gcd(A[*]), g[i] は g[i+1] を割り切る
例で見ましょう．
|2 6| → |2 6| → |2 2| → |2 0| |4 8| |0 -4| |0 4| |0 4|
つまり，上の A = [2,6;4,8] という行列は S A T = diag(2,4) となる，ということです．方程式 x A = b を解きたければ A を単因子標準系にして x S^{-1} = y とおいて y (S A T) = b T とし，y(i) = (B T)(i) / g(i) が解なのでこれを x = S y で戻せば良い，ということになります．これで一般の整数係数 A x = b が解けるようになりました．
実装では，基本変形は絶対値最小の要素を探してそれで行と列を掃き出すこと有限回で終わります（max |a*| が単調減少する）．これは１列の場合は Blankinship の手続きと一致するので，結局 Blankinship の拡張互除法は単因子標準形を求める掃き出し法の特殊版であったと理解できます．
一般に，ここまで述べたことはすべてユークリッド整域（ユークリッド互除法が定義できる領域）で成り立ちます．ユークリッド整域の例には他にも Z，Z[√-1], Z[√-2], Z[ω], K[x] などがあります．複素整数環 Z[√-1] や多項式環 K[x] の gcd は一度は書いておきたいかもしれません．
参考文献