拡張ユークリッド互除法によるベズーの等式の解の探索

はじめに

これは RSA完全理解 Advent Calendar の8日目の記事です。

さて、同カレンダーの過去の記事にて

ユークリッド互除法
ベズーの等式

について紹介したので、今回はユークリッド互除法をつかってベズーの等式の解を探索する方法について紹介します。
なお、このことを一般に拡張ユークリッド互除法などと呼んだりします。

拡張ユークリッド互除法とは

ユークリッド互除法は

ある値 $a, b$ についての割り算を $a = bq + r$ と表したときに
$gcd(a, b) = gcd(b, r)$ となる

という性質を利用して効率よく gcd(a, b) を求めるアルゴリズムのことでした。

前回の説明では素通りしましたが、じつは計算して得られる等式を整理すると、そのステップの余りを前のステップの情報をつかって表現できます！
以下はユークリッド互除法の手順と、そのステップの余りについて整理した式です。

$a\div b = q_1 \ ... \ r_1$ , 変形して $a = bq + r_1$ , さらに余りについて $a - bq = r_1$
$b\div r_1 = q_2 \ ... \ r_2$ , 変形して $b = rq_2 + r_2$ , さらに余りについて $b - rq_2 = r_2$
$r\div r_2 = q_3 \ ... \ r_3$ , 変形して $r = r_2q_3 + r_3$ , さらに余りについて $r - r_2q_3 = r_3$
以下割り切れるまで続ける...

となり、それぞれの割られる数を前のステップの値で埋めて、余りを表現することができます。

詳細はのちほど詳しく説明しますが、このことを利用すると最終的に $ax + by = gcd(a, b)$ の形になり、解である $x, y$ も一緒にわかるようになります。

当然 $gcd(a, b)$ は自身に1をかけた倍数であるので、

ベズーの等式:
$a, b$ を0でない整数としたとき、 $ax + by = d$ が解を持つ
$\iff$ $d$ は $gcd(a, b)$ の倍数

を満たす解のうちの一つです。

以上のように、ユークリッドの互除法の操作を応用すればベズーの等式の解を効率的に得る事ができます。

おお〜なるほどやってみようや

具体例あると理解が進むと思うので例をだしてやってみるぞ！

$a = 1152$ , $b = 504$ として普通にユークリッド互除法やってみる

$1152\div 504 = 2 \ ... \ 144$ , あまりについて整理して $144 = 1152 - 504 * 2$
$504\div 144 = 3 \ ... \ 72$ , あまりについて整理して $72 = 504 - 144 * 3$
$144\div 72 = 2 \ ... \ 0$ , あまりについて整理して $0 = 144 - 72 * 2$

となる。
$gcd(a, b)$ は割り切れたステップのひとつ前の余りであるため、そのときの $gcd(a, b)$ についての式

$72 = 504 - 144 * 3$

を出発点として前ステップの値で展開していく

$72 = 504 - 144 * 3$
$72 = 504 - (1152 - 504 * 2) * 3$ (144について前ステップの値を代入)
$72 = 504 - 3(1152) + 6(504)$ (↑を整理)
$72 = - 3(1152) + 7(504)$ (↑重複整理)

これで $-3\cdot 1152 + 7\cdot 504 = 72$ の式が得られ、1152と504の係数が $ax + by = gcd(a, b)$ の $x, y$ の解の一つとなります
念の為等式が成り立つか計算してみても

$72 = -3452 + 3528$

となり成立しています。

詳しい説明

まず、余り $r$ についての式をいくつか並べてみて、規則を見つけて一般化を試みます

$r_1 = a - bq_1$
$r_2 = b - rq_2$
$r_3 = r_1 - r_2q_3$
$r_4 = r_2 - r_3q_4$
$r_5 = r_3 - r_4q_5$
$...$

こうしてながめると 1, 2 の式以外は

$r_i = r_{i-2} - r_{i-1} q_i$

というような法則があることがわかります(ユークリッド互除法は商とあまりを引き継いでステップを進めるので、当然っちゃ当然)(ちなみに1, 2の式についてはうまく同じルールで取り扱うための手法をおまけで紹介してます)

さて、この式はじつは $r_i = ax_i + by_i$ の形で表すことができます

$r_i = r_{i-2} - r_{i-1} q_i$ は $r_i = ax_i + by_i$ とあらわせる

このとき $r_i = r_{i-2} + (- r_{i-1} q_i)$ と考えられるため、

$r_{i-1} = ma + nb$
$r_{i-2} = m'a + n'b$

であることが言えれば

$r_i = m'a + n'b + (- (ma + nb)q_i)$ $= a(m' - mq_i) + b(n' - nq_i)$

となり $r_i$ は $a, b$ の倍数の和であらわせます。

$r_i = r_{i-2} - r_{i-1} q_i$ から $r_i$ を表現するには直前の2つのステップの結果が必要であり、はじめの2つの式に $a, b$ の倍数の和である確認ができれば再帰的に適用することで後続も成り立つとわかります。

計算すると

$r_1 = a - bq_1 = a(1) - b(q_1)$ なので $r_1$ は $a, b$ の倍数の和
$r_2 = b - rq_2$
$= b - (a - bq_1)q_2$
$= b - aq_2 + bq_1q_2$
$= a(-1q_2) + b(1 + q_1q_2)$ なので $r_2$ は $a, b$ の倍数の和

となるので、後続の任意の $r_i$ について $a, b$ の倍数の和であることがわかります。
$a, b$ についての係数をそれぞれ $x_i, y_i$ とすると

$r_i = ax_i + by_i$

と表せます

これを元の式である $r_i = r_{i-2} - r_{i-1} q_i$ の $r_{i-2}$ と $r_{i-1}$ の部分に代入すると以下のようになります

$r_i = (ax_{i-2} + by_{i-2}) - (ax_{i-1} + by_{i-1})q_i \\ = ax_{i-2} + by_{i-2} - ax_{i-1}q_i - by_{i-1}q_i \\ = a(x_{i-2} - x_{i-1}q_i) + b(y_{i-2} - y_{i-1}q_i) \\ = ax_i + by_i$

これより $a, b$ の係数である $x_i, y_i$ について

$x_i = x_{i-2} - x_{i-1}q_i$
$y_i = y_{i-2} - y_{i-1}q_i$

といえることがわかります。
つまりユークリッド互除法において、手前2つのステップの $x, y$ の値と今のステップの商があれば $x_i, y_i$ が求められます

これをユークリッド互除法の処理が終了するまで、つまり $r_i = 0$ となるまで続けると、 $r_{i-1} = ax_{i-1} + by_{i-1} = gcd(a, b)$ となっているはずで、そのときの $a, b$ の係数である $x_{i-1}, y_{i-1}$ はベズーの等式 $ax + by = d$ を満たす解の一つとなっています！

おまけ: はじめの2ステップで参照される $a, b$ の取り扱いについて

ところで、ユークリッド互除法の最初の2ステップの処理は

$r_1 = a - bq_1$
$r_2 = b - rq_2$

のようになっておりそれぞれ $a, b$ を参照しています。

あまりが $a$ や $b$ となるステップは存在しないので、これらを参照する際には $r_i = r_{i-2} - r_{i-1} q_i$ の式における

$r_{i-2}$
$r_{i-1}$

が存在しません。

そこで $r_i = ax_i + by_i$ を守る形で $a$ と $b$ を表現するために、初期状態として以下の2項目を追加しておきます

$r_{-1} = a(1) + b(0)$ ( $x_{-1} = 1, y_{-1} = 0$ )
$r_0 = a(0) + b(1)$ ( $x_0 = 0, y_0 = 1$ )

このようにすると、はじめの2つのステップで $a, b$ について参照できます

これでユークリッド互除法を用いてベズーの等式 $ax + by = d$ を満たす解 $x, y$ が探索できることがわかりました！

実装

あとは説明したことを実装するだけです！

拡張ユークリッド互除法は以下のようにかけます(わかりやすさを優先して書き方は少し冗長です。ほんとうはもうちょっとシュッとかけます)

func extGCD(a, b int64) (gcd, x, y int64) {
    // はじめの2つのステップで a, b を参照するために初期状態を定義
    var prevx, prevy, curx, cury int64 = 1, 0, 0, 1
    for b != 0 {
        var q int64 = a / b
        a, b = b, a%b

        // x[i] = x[i-2] - q*x[i-1] から次のxを計算して prev, cur を進める
        nextx := prevx - q*curx
        prevx = curx
        curx = nextx

        // y[i] = y[i-2] - q*y[i-1] から次のyを計算して prev, cur を進める
        nexty := prevy - q*cury
        prevy = cury
        cury = nexty
    }
    return a, prevx, prevy
}