はじめに

導入

「ユークリッドの互除法」というのは、2つの整数の最大公約数(GCD)を計算する有名なアルゴリズムです。そして、これには「拡張ユークリッドの互除法」という応用があります。
今までこちらを何度も使ったことがあったのですが、あんまり深く気にしないで使っていたため、「どのような計算内容で」「なぜ成立するのか」を確かめて簡単にまとめました。

まとめ

先に今回の話のポイントです。

ユークリッドの互除法の応用に「拡張ユークリッドの互除法」がある。
拡張ユークリッドの互除法により、$ax+by=g~~( g=GCD(a,b) )$ の整数解 $(x,y)$ を求めることができる。
この「整数解を求めることができる」という性質により、互いに素な $(a,n)$ に対して、モジュロ逆数 $x\equiv a^{-1}~~mod~n$ の計算に使える。

ユークリッドの互除法のおさらい

ユークリッドの互除法は、2つの整数の最大公約数(GCD)を計算する高速なアルゴリズムです。
2数の割り算で余りを延々と計算していくことで、最終的には最大公約数まで辿り着くというものです。

例えば、$1095=3\cdot 5\cdot 73,~~474=2\cdot 3\cdot 79$ の2数の最大公約数は 3 です。
この2数に対し、以下のように得られた余りに差し替えて割り算を続けることで、最終的に割り切れるところまで行きますが、その直前の余りが最大公約数として得られるという寸法です。
※前回の割り算の除数 ( $\div$の右側 ) を次の被除数 ( $\div$の左側 )、余りを次の除数にする。

$$
\begin{array} \\
\color{purple}{1095}&\div&\color{blue}{474}&=&2&\dots&\color{red}{147} \\
\color{blue}{474}&\div&\color{red}{147}&=&3&\dots&\color{purple}{33} \\
\color{red}{147}&\div&\color{purple}{33}&=&4&\dots&\color{blue}{15} \\
\color{purple}{33}&\div&\color{blue}{15}&=&2&\dots&\underline{\color{red}{3}} ~\leftarrow GCD \\
\color{blue}{15}&\div&\color{red}{3}&=&5&\dots&0 \leftarrow End\\
\end{array}
$$

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

概要

では、本題の拡張ユークリッドの互除法です。
上でおさらいしたユークリッドの互除法では、2つの整数$(a,b)$の最大公約数$g$を求めることができました。
じつは、この最大公約数$g$に対して、2変数1次方程式 $ax+by=g$ は必ず$(x,y)$の整数解を持つことが分かっています。
例えば、上のおさらいで出した $(a,b)=(1095,474)~~( g=3 )$ であれば、$(x,y)=(-67,29),(91,-336),\cdots$ と、複数の解があります。
しかし、代表の解 $(x,y)=(-67,29)$ が求まれば、一般の解は $(x,y)=(-67+\frac{b}{g}\cdot k,~ 29-\frac{a}{g}\cdot k)$ ( $k$は任意の整数 ) としてまとめることができますので、1組解が見つかれば十分です。
元々のユークリッドの互除法に少し計算を付け加えることでこの「1組の解」を求められるのが、拡張ユークリッド互除法となっています。

不変条件の活用

では、拡張ユークリッド互除法の具体的な手順に進む前に、鍵となる部分を説明します。
それは、不変条件を活用することによる等式の生成です。

最終的な目標は $ax+by=g$ ですが、先に類似の等式 $ax+by=z$ という形を考えます。

そうしておいて、今次の2つの等式が成立しているものとします。

$ax''+by''=z''$
$ax'+by'=z'$

すると、なんらかの $q$ を持ち出すことで、次の等式を作り出せることが分かります。

$a(x''-qx')+b(y''-qy')=(z''-qz')$

これも $ax+by=z$ の形の等式であり、$(x,y,z)=(x''-qx',y''-qy',z''-qz')$ のケースに該当するものです。

であれば、$ax+by=z$ の形の等式を次々作り出し、$z=g$ になるように $z$ を変化させることができれば、最終目標を達成することができます。これが拡張ユークリッドの互除法の鍵となっています。

アルゴリズムの構成

では、具体的なアルゴリズムの構成に移ります。
上で「$z=g$ になるように $z$ を変化させることができれば」という話をしました。ここに、元々のユークリッドの互除法を併せてみます。

ユークリッドの互除法は、言うなれば次のように $z$ を変化させる ( 数列 ${z_k}$ を構成する ) 計算法でした。

$z_0=a$ とする
$z_1=b$ とする
$k\ge 2$ に対し、$z_k=mod(z_{k-2},z_{k-1})$ ( $z_{k-2}\div z_{k-1}$ の余り ) とする

この最後を次のように言い換えます。

$k\ge 2$ に対し、$q_k=\lfloor \frac{z_{k-2}}{z_{k-1}} \rfloor$ とする
※この $q_k$ は $z_{k-2}\div z_{k-1}$ の商にあたる数値です
$k\ge 2$ に対し、$z_k=z_{k-2}-q_{k}z_{k-1}$ とする

これは、余りの計算を商を使って言い換えたものです。
こうすることで、前章で挙げたような $(z'',z')\rightarrow z''-qz'$ の変形に合わせることができていることが分かります。
これに、さらに $x,y$ を併せます。

$x_0=1, y_0=0, z_0=a$ とする。
この時 $ax_0+by_0=z_0$ が成立する。
$x_1=0, y_1=1, z_1=b$ とする。
この時 $ax_1+by_1=z_1$ が成立する。
$k\ge 2$ に対し、$q_k=\lfloor \frac{z_{k-2}}{z_{k-1}} \rfloor$ とする
$k\ge 2$ に対し、$x_k=x_{k-2}-q_kx_{k-1},~~y_k=y_{k-2}-q_ky_{k-1},~~z_k=z_{k-2}-q_kz_{k-1}$ とする。
この時、$x_{k-1},y_{k-1},z_{k-1}$ までと同様に $ax_k+by_k=z_k$ が成立する。
この ${z_k}$ はユークリッドの互除法の計算の結果をなぞる数列となっており、いつか $g$ に到達する。
その時 ( $z_k=g$ の時 ) の $(x_k,y_k)$ が、$ax+by=g$ の解 $(x,y)$ の1組である。

これが拡張ユークリッドの互除法の計算内容です。

実施例

それでは、実際に $(a,b)=(1095,474)~~(g=3)$ に拡張ユークリッドの互除法を適用した時の、各$k$ に対する $x_k,y_k,z_k,q_k$ の値の変化を実際に見てみると、次の表のようになります。

$k$	$z_k$	$q_k=\lfloor \frac{z_{k-2}}{z_{k-1}}\rfloor$	$x_k$	$y_k$	$ax_k+by_k$
0	1095	-	1	0	1095
1	474	-	0	1	474
2	147	2	1	-2	147
3	33	3	-3	7	33
4	15	4	13	-30	15
5	3	2	-29	67	3
6	0	5	-	-	-

例えば $k=3$ のところですが、以下のような計算になっています。

2つ前の項 $(x_1,y_1,z_1)=(0,1,474)$
1つ前の項 $(x_2,y_2,z_2)=(1,-2,147)$
$z_k$ の割り算の結果として $q_3=\lfloor \frac{474}{147} \rfloor=3$
$x_3=x_1-3x_2=-3,~~y_3=y_1-3y_2=7,~~z_3=z_1-3z_2=33$

そうして、いずれの $k$ でも $ax_k+by_k=z_k$ が成立しており、$z_k=g~(=3)$ となる $k=5$ の $(x_5,y_5)=(-29,67)$ これが $ax+by=g$ の解として求められている、ということになります。

なお、この計算方法のRubyでの実装例は次の通りです。
※簡単のため、正整数のみ受け付ける実装としています。

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

def extended_eucledean(a,b)
  raise unless a>0 && b>0
  xp,xc = 1,0
  yp,yc = 0,1
  zp,zc = a,b
  loop do
    q = zp/zc
    zp,zc = zc,zp-q*zc
    break if zc==0
    xp,xc = xc,xp-q*xc
    yp,yc = yc,yp-q*yc
  end
  return xc,yc,zp
end

上記実装で $(a,b)=(1095,474)$ のケースを試すと、ちゃんと上に挙げた計算結果が得られ、それが正しいことが分かります。

実行例

irb(main):002:0> a,b = 1095,474
=> [1095, 474]
irb(main):003:0> x,y,g = extended_eucledean(a,b)
=> [-29, 67, 3]
irb(main):004:0> a.gcd(b) == g and puts "GCD(#{a},#{b}) = #{g}"
GCD(1095,474) = 3
=> nil
irb(main):005:0> a*x+b*y == g and puts "#{a}*(#{x})+#{b}*(#{y})=#{g}"
1095*(-29)+474*(67)=3
=> nil

応用

モジュラ逆数計算

最後におまけてきになりますが、そもそも「なぜ拡張ユークリッドの互除法の話を整理したか」という動機に関わる話です。
それは、( 一般の話としてどうかは分かりませんが ) 「モジュラ逆数計算」に活用できるからです。

モジュラ逆数は、何らかの$n$を法とするmod計算の中で、何らかの数$a$に対して、$ax\equiv 1~mod~n$ が成立する $x$ のこと、$a^{-1}$ を指します。
すなわち、mod計算の中での逆数であり、割り算に相当する計算の鍵となるものです。
※$a,n$の条件によっては、モジュラ逆数が存在しない場合もあります。

modでない計算であれば、例えば $3x=1$ の逆計算、割り算で $x=3^{-1}=\frac{1}{3}$ と計算できます。
しかし、$3x\equiv 1~mod~7$ のような mod計算では同じような計算は通じません。
※この例では $3\cdot 5\equiv 1~mod~7$ のため、$3^{-1}\equiv 5~mod~7$ となります。

このモジュラ逆数計算は、私がよく題材として取り上げる公開鍵暗号でも、RSA暗号/署名での $ed\equiv 1~mod~(p-1)(q-1)$ から $d$ を求める計算や、ECDSA(署名)の内部での計算等、様々な場面で現れます。そのため、効率的な計算方法は非常に重要なのです。

拡張ユークリッドの互除法のモジュラ逆数計算への応用は単純です。
2整数 $a,n$ が互いに素 ( 最大公約数が 1 ) なのを前提として、この組み合わせに拡張ユークリッドの互除法を適用します。
そうすると、得られる $(x,y)$ は $ax+ny=1 \Rightarrow ax\equiv 1~mod~n$ を満たすため、$x$ の方がそのままモジュラ逆数 $a^{-1}$ になっているという寸法です。

おわりに

有名なアルゴリズムなので、今更な話に感じる人もいるかもしれませんが、アルゴリズムの裏に隠れている仕組みの把握に役立てば幸いです。