[数学] Euclid の互除法

$x,y$ の最大公約数 $\gcd (x,y)$ を計算します（最大公約数は英語でgreatest common divisorですので、gcdと略されます）。

なお、$x=y=0$ の場合には、任意の整数が $x,y$ の公約数となります。そのため公約数の最大値は存在せず、$\gcd (0,0)$は通常定義しません（イデアル論的な観点から、$0$と定義することもあります。）

数学的背景

$d$ が $n$ を割り切るとき、$d\mid n$ と書きます。次が成り立ちます：
$$d\mid x\text{かつ}d\mid y⟺d\mid x\text{かつ}d\mid x–ky$$

これは $d$ の倍数同士の和・差が $d$ の倍数であることから、簡単に証明できます。さらにこの定理から、$\gcd$ の性質$$\gcd (x,y)=\gcd (x,y-kx)$$が成り立つことが分かります。このことを用いると、$2$ 数の最大公約数を次のように計算できます：
$$\gcd (345,100)=\gcd (100,45)=\gcd (45,10)=\gcd (10,5)=\gcd (5,0)=5.$$ この計算方法を、Euclidの互除法といいます。大きい数を、小さい数で割った剰余に置き換えることを繰り返しています。

計算量

$a>b>0$ であり，$a$ を $b$ で割った余りが $c$ であるとします。
$$a\geqq 2b⟹b+c\leqq b+b\leqq \frac{2}{3}(a+b),a\leqq 2b⟹b+c\leqq b+(a-b)\leqq \frac{2}{3}(a+b)$$より、Euclidの互除法の計算を1ステップ進めると、2数の和が $\frac{2}{3}$ 倍以下になることが分かります。したがって、$\gcd (x,y)$ の計算は $O(\log (x+y))$ で行えます。1回目の除算でどちらも $min(x,y)$ 以下になるため、$O(\log min(x,y))$ とすることも可能です。（上の証明から、計算回数の上界が ${\log }_{1.5}(x+y)$ と評価できますが，この底 $1.5$ は精密化すると，黄金比$\frac{1+\sqrt{5}}{2}=1.618$ まで改善できます。）

我々の見慣れている $2,3$ 桁の整数であれば、素因数から最大公約数を求めるのが簡潔そうに見えますが、桁数が大きくなるとEuclidの互除法の方が圧倒的に性能が高く、計算機では数百桁といった数のgcdですら瞬時に計算することができます。

実装

gcd関数は、math・numpyなどのライブラリにも備わっていますが、ここでは自作してみます。

def gcd(a, b):
    while b:
        a, b = b, a % b
    return a

def gcd(a, b):
    return gcd(b, a % b) if b else a

ループによる実装、再帰による実装の2つを挙げていますが、どちらも処理内容は同じです。
なお、関数呼び出しのオーバーヘッドの分、ループによる実装の方がパフォーマンスは良いです。
（関数の1回あたりの処理がある程度重いならば、呼び出しの負荷は相対的に無視量になりますが、今回は1回除算をする程度ですので、関数の呼び出しの差が相対的に多くなります。パフォーマンスについては別記事で、計測結果などを扱っていく予定です。）

関連：最小公倍数

$2$ 数 $a,b$ の最小公倍数を、$\mathrm{lcm}(a,b)$ と書きます（least common multiple）。$\gcd$ との間に次の関係があります：
$$\mathrm{lcm}(a,b)=\frac{ab}{\gcd (a,b)}.$$ このことは、素因数分解を考えることで証明されます（3数以上の最大公約数・最小公倍数の間にはこのような関係式はありません）。

したがって、$\mathrm{lcm}(a,b)$ も $O(\log (a+b))$ で計算できることが分かります。

def lcm(a,b):
    return a // gcd(a,b) * b

関連：ax+by=dの整数解

$d=\gcd (a,b)$ とするとき、$ax+by=d$ の解の1つをEuclidの互除法により計算できます。このことは高校数学でもしばしば取り上げられます。「拡張互除法」ということもあるようです。

def extgcd(a, b):
    """
    ax + by = gcd(a,b) = d となる (x,y,d) を返す
    """
    if b == 0:
        return (1, 0, a)
    q, r = a // b, a % b
    x, y, d = extgcd(b, r)
    s, t = y, x - q * y
    return s, t, d  

関連：a mod Nの逆元

$a,N$ が整数で、$\gcd (a,N)=1$ のとき、$ax\equiv 1(\mathrm{mod}N)$ となる $x$ を $amodN$ の逆元といいます。上で述べたように、Euclidの互除法から $ax+Ny=1$ を満たす整数 $(x,y)$ が計算できます。この $x$ は当然 $ax\equiv 1(\mathrm{mod}N)$ を満たすので、これが逆元です。$x<0$ になって欲しくない場合は、負数ならば $N$ を加える処理を追加すればよいです。

def inv_mod(a, N):
    """
    ax = 1 mod N となる x を返す
    """
    return extgcd(a, N)[0]

関連：中国剰余定理

$a_1$, $a_2$を整数として、$m_1$, $m_2$ を互いに素な整数とします。このとき
$$x\equiv a_1(\mathrm{mod}{m}_1),x\equiv a_2(\mathrm{mod}{m}_2)$$を満たす整数 $x$ が、$mod{m}_1{m}_2$ で一意に定まります。このことを中国剰余定理といいます。

$x$ を計算してみましょう。第1の条件が成り立つような $x$ は、$x=a_1+m_{1}y$ と書けます。上手く $y$ を調整して、$x\equiv a_2(\mathrm{mod}{m}_2)$ が成り立つようにしましょう。
$$x\equiv a_2(\mathrm{mod}{m}_2)⟺a_1+m_{1}y\equiv a_2(\mathrm{mod}{m}_2)⟺m_{1}y\equiv a_2-a_1(\mathrm{mod}{m}_2)$$ となります。したがって、$m_{1}mod{m}_2$ の逆元を $z$ とすれば、$y=(a_2–a_1)z$ が条件を満たします。

def CRT(a_1, a_2, m_1, m_2):
    """
    x = a_1 mod m_1, x = a_2 mod m_2 となる x を返す
    """
    y = (a_2 - a_1) * inv_mod(m_1, m_2) % m_2
    return a_1 + m_1 * y