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

# はじめに
## 導入
「ユークリッドの互除法」というのは、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での実装例は次の通りです。
※簡単のため、正整数のみ受け付ける実装としています。

```rb:拡張ユークリッドの互除法実装
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:実行例
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}$ になっているという寸法です。

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