整数論１ ユークリッド互除法

　ユークリッド互除法は２つの整数の最大公約数を求めるアルゴリズムである．証明も含めて紹介する．

定理１： を整数，をで割った余りをとする．このとき$$\mathrm{gcd}(a,b) = \mathrm{gcd}(b, r)$$である．

証明：

とする．定義よりを満たす整数が存在する．従ってはの公約数であり．またよりはの公約数である．従って．これよりとなる．（証明終了）

この定理１を下にユークリッド互除法を説明する．2つの自然数の最大公約数を求めたい．数列を次のように定める

1. としをで割った余りをとする．
2. をで割った余りをとする．

ここではになるまで単調減少である．で初めてになるとすると定理1より$$\mathrm{gcd}(a,b) =\mathrm{gcd}(r_0,r_1)= \dots \mathrm{gcd}(r_{n_0 - 1},r_{n_0})=\mathrm{gcd}(r_{n_0 - 1}, 0) = r_{n_0 - 1}$$となる．これにより最大公約数を求めることができる．これがユークリッド互除法である．

　ユークリッド互除法の過程を詳しく見ることで次の定理が従う．

定理２ をでない整数とする．このとき$$ax+by=\mathrm{gcd} (a,b)$$を満たす整数が存在する．この等式をべズーの等式という．

証明： と仮定してよい．実際に対しを求め，符号を調節すればよい．としてを定める．このとき \begin{equation*} \begin{pmatrix} r_0\\ r_1 \end{pmatrix} = \begin{pmatrix} k_0 & 1 \\ 1 & 0 \end{pmatrix} \cdots \begin{pmatrix} k_{n_{0}-1} & 1\\ 1 & 0 \end{pmatrix} \begin{pmatrix} r_{n_0} \\ 0 \end{pmatrix} \end{equation*} とかける．行列$$\begin{pmatrix}k_n & 1\\ 1 & 0 \end{pmatrix}$$は行列式より逆行列を持つ．従ってについて解くことができる．より定理が示された．（証明終了）

# coding: utf-8
import numpy as np
# ユークリッド互除法

def Euclid(a, b):
    # ユークリッド互除法を用いて最大公約数を求める関数
    a = abs(a)
    b = abs(b)
    if a == 0 and b == 0: raise ValueError("両方0です")
    elif a == 0: return b
    elif b == 0: return a
    else:
        a, b = max(a, b), min(a, b)
        while b != 0:
            a, b = b, a % b
        return a

def Bezout_identity(a, b):
    # べズーの等式 ax + by = gcd(a,b) を満たすx,yを求める関数

    if a == 0 and b == 0: raise ValueError("両方0です")
    elif a == 0: return 0, 1 
    elif b == 0: return 1, 0
    else:
        c = a
        d = b
        k = np.array([[1, 0], [0, 1]])  # 単位行列
        a = abs(a)
        b = abs(b)
        junban = a > b                  # a,bの大小関係を記憶
        a, b = max(a, b), min(a, b)
        while b != 0:
            v = np.array([[a // b, 1], [1, 0]])
            k = np.dot(k,v)
            a, b = b, a % b
        k = np.linalg.inv(k)
        x = int((c // abs(c)) * round(k[0][0]))
        y = int((d // abs(d)) * round(k[0][1]))
        if junban == True:
            ans = c * x + d * y
            if ans > 0: return x, y
            else      : return -x, -y
        else:
            ans = c * y + d * x
            if ans > 0: return y, x
            else      : return -y, -x
        
if __name__ == "__main__":
    for (a, b) in [(47, 18), (-14, 26)]:
        print("###########################")
        print(a, b)
        print("最大公約数{0}".format(Euclid(a, b)))
        print(Bezout_identity(a, b))

###########################

47 18

最大公約数1

(5, -13)

###########################

-14 26

最大公約数2

(-2, -1)