拡張ユークリッド互除法の話

先日の連分数の計算と拡張ユークリッド互除法って少し似てるなと思ったというメモ。（連分数の計算で次の項を求める時にやっている操作ってユークリッドの互除法そのものだし当然か）

拡張ユークリッド互除法というのは、 $ax+by=gcd(a,b)$ を満たす $x,y\in\mathbb{Z}$ を求める方法。
$a,b$ を $gcd(a,b)$ で割れば $a'x+b'y=1$ （a'とb'は互いに素）という形になる。

（以下、 $a\lt b$ の場合 swap(a, b)してから）
ユークリッドの互除法の各ステップが
$\begin{pmatrix} a_0 \\ b_0 \end{pmatrix}=\begin{pmatrix} a \\ b \end{pmatrix}, \begin{pmatrix} q_0 \\ r_0 \end{pmatrix}=\begin{pmatrix} a_0\backslash b_0 \\ a_0\%b_0 \end{pmatrix},$
$\begin{pmatrix} a_1 \\ b_1 \end{pmatrix}=\begin{pmatrix} b_0 \\ r_0 \end{pmatrix}, \begin{pmatrix} q_1 \\ r_1 \end{pmatrix}=\begin{pmatrix} a_1\backslash b_1 \\ a_1\%b_1 \end{pmatrix},$
$\cdots,$
$\begin{pmatrix} a_{k-1} \\ b_{k-1} \end{pmatrix}=\begin{pmatrix} b_{k-2} \\ r_{k-2} \end{pmatrix}, \begin{pmatrix} q_{k-1} \\ r_{k-1} \end{pmatrix}=\begin{pmatrix} a_{k-1}\backslash b_{k-1} \\ a_{k-1}\%b_{k-1} \end{pmatrix},$
$\begin{pmatrix} a_k \\ b_k \end{pmatrix}=\begin{pmatrix} b_{k-1} \\ r_{k-1} \end{pmatrix}=\begin{pmatrix} gcd(a,b) \\ 0 \end{pmatrix}$
のとき
$\begin{pmatrix} x \\ y \end{pmatrix}= \begin{pmatrix} 0&1 \\ 1&-q_0 \end{pmatrix}\begin{pmatrix} 0&1 \\ 1&-q_1 \end{pmatrix} \cdots\begin{pmatrix} 0&1 \\ 1&-q_{k-1} \end{pmatrix}\begin{pmatrix} 1 \\ 0 \end{pmatrix}$
という風に $x,y$ が求まる、というのが拡張ユークリッド互除法である。

ちなみに $\displaystyle\frac{a}{b}$ は連分数で $\displaystyle[q_0;q_1,q_2,\cdots,q_{k-1}]=q_0+\frac{1}{q_1+\frac{1}{q_2+\frac{1}{\cdots}}}$ と表すことができる。

// extgcd.cc
#include <cstdio>
#include <cstdlib>

int extgcd(int a, int b, int &x, int &y, int c=1, int d=0, int e=0, int f=1) {
  if (b == 0) {
    x = c; y = e;  // dot([c,d; e,f], [1;0])
    return a;
  }
  int q = a/b, r = a%b;
  return extgcd(b, r, x, y, d, c-q*d, f, e-q*f);  // dot([c,d; e,f], [0,1; 1,-q])
}

int main(int argc, char **argv) {
  if (argc == 3) {
    int a = atoi(argv[1]), b = atoi(argv[2]);
    int x, y;
    int g = extgcd(a, b, x, y);
    printf("gcd(%d,%d)=%d; x=%d,y=%d\n", a, b, g, x, y);
  } else {
    printf("usage: %s <a> <b>\n", argv[0]);
  }
  return 0;
}

おまけ

in scheme

(define (extgcd a b)  ; -> gcd(a,b), x, y
  (let loop ((a a) (b b) (c 1) (d 0) (e 0) (f 1))
    (if (zero? b) (values a c e)
        (let ((q (quotient a b)))
          (loop b (remainder a b) d (- c (* q d)) f (- e (* q f)))))))

(receive (g x y) (extgcd 4 3)
  (format #t "gcd(a,b)=~d; x=~d,y=~d\n" g x y))
;; => gcd(a,b)=1; x=1,y=-1

(receive (g x y) (extgcd 48 15)
  (format #t "gcd(a,b)=~d; x=~d,y=~d\n" g x y))
;; => gcd(a,b)=3; x=1,y=-3

↑商と剰余の計算をGaucheだと

(receive (q r) (quotient&remainder a b)
  (loop b r d (- c (* q d)) f (- e (* q f))))

みたいに書けるけど他の処理系でも行ける？

naoya_t@hatenablog

いわゆるチラシノウラであります

拡張ユークリッド互除法の話

おまけ