feature image

2021年4月18日 | ブログ記事

ベズー係数とN項の拡張ユークリッドの互除法

こんにちは、すく(0214sh7)です
東工大数学系の3年生で、traPではアルゴリズム班のリーダーをやっています。

本来この時期にブログを書くつもりはなかったのですが、唐突にベズー係数すきすき侍と化したので筆をとりました。

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

競技プログラミングをやっている中で拡張ユークリッドの互除法(extended Euclidean algorithm)の存在を知った人は多いと思います。
機能としては通常のユークリッドの互除法の拡張になっています。
けんちょんさんのqiita記事などで丁寧に解説されていますので、この記事では詳しく解説しません。

ベズー係数

これはベズーの等式(Bézout's identity)と呼ばれる式です。
フランスの数学者エティエンヌ・ベズーからこの名前がつきました。
ベズーの等式にまつわる定理として、ベズーの補題(Bézout's lemma)があります。この定理は以下を主張しています

を0でない整数、とする。このとき、1次方程式

は、の倍数であれば無数の整数解を持ち、の倍数でなければ整数解を持たない。

証明(クリックで開く)

まずはの倍数でなければ整数解は存在しないことを示したい。
の倍数であるので、整数が存在してである。
同様に、整数が存在してである。
すると、である。が整数だとすればも整数であるので、の倍数である。
このため、の倍数でなければを満たす整数は存在しない。

次にの倍数であれば整数解を1つ持つことを示したい。
ここで、を満たす整数が存在することを用いる。補題として以下に示す。

補題の証明(クリックで開く)

一般性を失わずとする。
として、をみたすである限り存在する。そこで、と定義する。
をみたす最小の(存在する)をとすると、である。
これより、をみたすとしてが存在することがわかる。

をみたすが存在すると仮定すると、より、をみたすとしてが存在する。
これを再帰的に適用することで、を満たすの存在が言える。


また、の倍数であるので、整数が存在してである。
すると、整数が存在しを満たすため、確かに整数解としてふさわしい。

そして整数解が1つあれば無数にあることを示したい。
これは難しくない。が整数解であればであるのでこれを変形して。よってが整数解になる。
これを好きなだけ繰り返すことでいくらでも異なる整数解を生成できる。よって整数解は無数にあると言える。

証明は省略するが、特殊解と整数を用いてと表せる解が、方程式の解のすべてである。


このことからに整数解が存在するの中で絶対値最小であり、任意のの倍数である、ということが言えます。

ということなので、をみたすベズー係数(Bézout coefficients)と名付けます(一意ではないことに気をつけてください)。

このベズー係数を1つ見つけることができるアルゴリズムとして知られているのが、拡張ユークリッドの互除法なのです!
最大公約数を求めるアルゴリズムであるユークリッドの互除法に機能を追加してベズー係数を求められるようにしていることから、この名前がつきました。

ここで驚くべき性質があります。拡張ユークリッドの互除法が与えるこのベズー係数はなんと、

をみたす2つのうち1つであるのです!

証明(クリックで開く)

帰納法を用いる。

一般性を失わずとする。
まず拡張ユークリッドの互除法を実行してに至る直前のステップ、つまりあるがあってという場合を考える。
この場合与えられるベズー係数はであるので、よりOK。

次にを満たすが条件を満たしていたと仮定した場合にを満たすが条件を満たすこと、すなわちとして

を示す。
より、

より、

以上より、機能法を用いて拡張ユークリッドの互除法が与えるベズー係数が

を満たすものであることが示された。


2項のベズー係数を拡張ユークリッドの互除法を用いて計算できること、わかっていただけましたか?
C++による実装例はこちらから見れます。

3項でのベズー係数

ここからが本題です
N項での議論をする前に前準備として3項での議論をしようと思います。
実はベズー係数は3項としてもその存在が言えます。つまり、を0でない整数として

をみたす整数(x,y,z)が存在します(そしてはそれを存在たらしめる右項のうち絶対値最小)。

これは次のようにして存在証明・計算ができます。

  1. のベズー係数を拡張ユークリッドの互除法を用いて一つ求める。
  2. のベズー係数を拡張ユークリッドの互除法を用いて一つ求める。
  3. のベズー係数の一つである。

これが何故かというと、
の定義より、が成り立ち、2つ目の式の最初の部分に1つ目の式を代入することで、

を導けるからです。
なので、はベズー係数の定義を満たしています。

N項でのベズー係数

書いてて思ったんだけどニッチすぎませんか?競プロですら出てくるか怪しいぞ

3項のものと同様の手順でN項に拡張できます。
最後の1項を除くN-1項のベズー係数が見つかっているものとします。
つまり、のベズー係数の一つであることがわかっているものとします。また、とします。

  1. のベズー係数を拡張ユークリッドの互除法を用いて一つ求める
  2. のベズー係数の一つである。

3項のものと同じ理屈でこれでベズー係数を求められることがわかります。
のときはを満たせばよいので、ベズー係数はとしておくとよいでしょう。

これをするコードを書いてみましょう。

ソースコード

import math

def bezout(a,b):
    if b==0:
        return [1,0]
    r = a%b
    q = (a-r)//b
    
    d = bezout(b,r)
    D = [d[1],d[0]-q*d[1]]
    return D
    
def bezoutN(A):
    if len(A)<=0:
        return []
    X = [1]
    g = A[0]
    for i in range(len(A)):
        if i==0:
            continue
        B = bezout(g,A[i])
        p = B[0];q = B[1]
        for j in range(len(X)):
            X[j] = p*X[j]
        X.append(q)
        g = math.gcd(g,A[i])
    return X


print(bezout(111,30))
print(bezoutN([111,30,7]))
print(bezoutN([111,30,7,11]))

出力

[3, -11]
[-6, 22, 1]
[-6, 22, 1, 0]

いい具合ですね。

bezoutNの時間計算量はです。
ところで、このコードにはいくつか改善できる点があります。

時間計算量改善

というと、パット見がネックになってそうです。
これはベズー係数を保存しておいて、掛け算を後でまとめてやるだけでOKです。
変更後のオーダーはです。

ソースコード

import math

def bezout(a,b):
    if b==0:
        return [1,0]
    r = a%b
    q = (a-r)//b
    
    d = bezout(b,r)
    D = [d[1],d[0]-q*d[1]]
    return D
    
def bezoutN(A):
    if len(A)<=0:
        return []
    g = A[0]
    P = [];Q=[]
    for i in range(len(A)):
        if i==0:
            continue
        B = bezout(g,A[i])
        p = B[0];q = B[1]
        P.append(p)
        Q.append(q)
        g = math.gcd(g,A[i])
    
    X = [0]*(len(A))
    now = 1
    for i in reversed(range(len(A)-1)):
        X[i+1]=now*Q[i]
        now = now*P[i]
    X[0]=now
    return X


print(bezout(111,30))
print(bezoutN([111,30,7]))
print(bezoutN([111,30,7,11]))

出力

[3, -11]
[-6, 22, 1]
[-6, 22, 1, 0]

オーバーフロー対策その1

この計算だと、項のうち前の方、とくに初項は個のベズー係数がかかっています。
オーバーフローを心配したくなりますよね。
そこで、前の係数をできるだけ小さくしたいです。

ここでベズーの等式の性質を使います。
として、

をみたすを拡張ユークリッドの互除法が与えてくれるということは先述しました。
なので、が最小であるようなベズー係数は拡張ユークリッドの互除法がもたらしたものをとして、

のどこかにあるということになります。
が最小なベズー係数を見つけてあげれば、オーバーフローを起こしにくくなります。(気休め程度な気もするが・・・)

オーバーフロー対策その2

気休め程度が気に食わないので、オーダーレベルで対策をします。

を合体すると、

を得ます。

つまり、のベズー係数のベズー係数がわかっていれば、

  1. のベズー係数を拡張ユークリッドの互除法を用いて一つ求める。
  2. のベズー係数の一つである。

という手順でのベズー係数がで計算できます。
ではをどうやって計算するかというと、再帰を使います。
セグ木やFFTと同じですね。単項になったらそのベズー係数はでしょう。

これをうまく使うことででN項のベズー係数を計算できます。logが一つ余計につきましたが、加算なのでそれもあまり気になりません。それに、各項にかかるベズー係数の個数が個になってます!これはかなりオーバーフロー対策と言えるのではないでしょうか?
実はレベルの係数を与えられると2,3発でアウトなのだが・・・

ソースコード

import math

def bezout(a,b):
    if b==0:
        return [1,0]
    r = a%b
    q = (a-r)//b
    
    d = bezout(b,r)
    D = [d[1],d[0]-q*d[1]]
    return D
    
def bezoutN(A):
    if len(A)<=0:
        return []
    if len(A)==1:
        return [1]
    A1 = A[:(len(A)//2)]
    A2 = A[(len(A)//2):]
    
    #実はここのgcdを可読性のために愚直計算にしている影響で計算量が
    #O(Nlog(max(A))+NlogN^2)になり、logが一つ増えます
    #bezoutNがgcdを同時に返すことでこれは解消します
    g1 = A1[0]
    for i in range(len(A1)):
        g1 = math.gcd(g1,A1[i])
    g2 = A2[0]
    for i in range(len(A2)):
        g2 = math.gcd(g2,A2[i])
    
    B = bezout(g1,g2)
    p = B[0];q=B[1]
    x = bezoutN(A1)
    y = bezoutN(A2)
    
    X = [0]*(len(A))
    for i in range(len(x)):
        X[i] = x[i]*p
    for i in range(len(y)):
        X[len(A1)+i] = y[i]*q
    
    return X


print(bezout(111,30))
print(bezoutN([111,30,7]))
print(bezoutN([111,30,7,11]))

出力

[3, -11]
[0, -3, 13]
[0, 0, -3, 2]

出力されているベズー係数が変わりましたね!
2項のときに一意でないと言ったように、N項でもベズー係数は一意ではありません。ですがからわかる通り、ベズー係数としての定義は満たしています!

おまけ:中国剰余定理

競技プログラマーの間では単に中国剰余定理あるいはCRTと呼ばれるテクニックである、
互いに素な整数と整数があり、

であるとき

を満たす整数を求める。
というアルゴリズムは、拡張ユークリッドの互除法を用いて実行できます!

  1. のベズー係数を拡張ユークリッドの互除法を用いて一つ求める
  2. (をで割ったあまり)である。

あとがき

最初はextGCDを自作ライブラリに追加するために拡張ユークリッドの互除法の勉強を始めたんですよ。

extGCD持ってなくて必要になった時いつもけんちょんのブログからパクってたんだけどとうとう自前のextGCDが完成した

— ⚡すく⚡ (@0214sh7) April 7, 2021

それで次はCRTの勉強・ライブラリ化に行くのかと思いきやライブラリページの命名に悩んだ結果ベズー係数の存在を知り、気づけばベズー係数すきすき侍と化していました。
数学科に入った結果代数への解像度が上がったせいです 数学科ってこわいね
このブログも代数の教科書を見ながら書いています。

そしてN項のベズー係数を求める方法を探した結果日本語の記事が見つからなかった(英語はあった)ので、この記事を作成しました。
N項のベズー係数を求める方法を知りたい人に届くことを祈っています。

参考文献

Harold M. Stark, An Introduction to Number Theory,MIT press,The MIT Press,1978.(芹沢正三,安藤四郎 訳『初等整数論』,現代数学社).

中島匠一,『代数と数論の基礎』,共立出版,2000.

ベズーの等式,Wikipedia,2020年11月18日更新,2021年4月14日アクセス.
https://ja.wikipedia.org/wiki/ベズーの等式

Extended Euclidean algorithm,Wikipedia,2021年1月29日更新,2021年4月14日アクセス.
https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm

けんちょん (Otsuki),拡張ユークリッドの互除法 〜 一次不定方程式 ax + by = c の解き方 〜,Qiita,2019年10月17日更新,2021年4月14日アクセス.
https://qiita.com/drken/items/b97ff231e43bce50199a

0214sh7 icon
この記事を書いた人
0214sh7

あおこーだー アルゴリズム班の班長さんです

この記事をシェア

このエントリーをはてなブックマークに追加
共有

関連する記事

2017年11月14日
IBIS2017参加報告
Keijan icon Keijan
ERC20トークンを用いた宝探しゲーム(真)の提案【アドベントカレンダー2018 10日目】 feature image
2018年11月3日
ERC20トークンを用いた宝探しゲーム(真)の提案【アドベントカレンダー2018 10日目】
Azon icon Azon
2018年12月12日
多重スリーブの世界と,各種進捗報告。
Silviase icon Silviase
2021年4月26日
CPCTF2021 作問者writeup by hukuda222
hukuda222 icon hukuda222
2021年4月25日
CPCTF2021 作問者writeup by xxpoxx
xxpoxx icon xxpoxx
2018年12月16日
ICPCアジア地区横浜大会参加記【アドベントカレンダー2018 52日目】
eiya icon eiya
記事一覧 タグ一覧 Google アナリティクスについて