こんにちは、すく(0214sh7)です
東工大数学系の3年生で、traPではアルゴリズム班のリーダーをやっています。
本来この時期にブログを書くつもりはなかったのですが、唐突にベズー係数すきすき侍と化したので筆をとりました。
拡張ユークリッドの互除法とは?
競技プログラミングをやっている中で拡張ユークリッドの互除法(extended Euclidean algorithm)の存在を知った人は多いと思います。
機能としては通常のユークリッドの互除法の拡張になっています。
けんちょんさんのqiita記事などで丁寧に解説されていますので、この記事では詳しく解説しません。
ベズー係数
これはベズーの等式(Bézout's identity)と呼ばれる式です。
フランスの数学者エティエンヌ・ベズーからこの名前がつきました。
ベズーの等式にまつわる定理として、ベズーの補題(Bézout's lemma)があります。この定理は以下を主張しています
を0でない整数、をとする。このとき、1次方程式
は、がの倍数であれば無数の整数解を持ち、がの倍数でなければ整数解を持たない。
証明(クリックで開く)
まずはがの倍数でなければ整数解は存在しないことを示したい。
はの倍数であるので、整数が存在してである。
同様に、整数が存在してである。
すると、である。が整数だとすればも整数であるので、はの倍数である。
このため、がの倍数でなければを満たす整数は存在しない。
次にがの倍数であれば整数解を1つ持つことを示したい。
ここで、を満たす整数が存在することを用いる。補題として以下に示す。
補題の証明(クリックで開く)
一般性を失わずa,b>0とする。
a0:=a,b0:=b、i≥0として、ai=biqi+ri(0≤ri<bi)をみたすqi,riがbi=0である限り存在する。そこで、ai+1:=bi,bi+1:=riと定義する。
bi=0をみたす最小のi(存在する)をkとすると、ak=dである。
これより、akx+bky=dをみたす(x,y)として(1,0)が存在することがわかる。
ai+1z+bi+1w=dをみたす(z,w)が存在すると仮定すると、ai+1z+bi+1w=biz+riw=biz+(ai−biqi)w=aiw+bi(z−qiw)=dより、aix+biy=dをみたす(x,y)として(w,z−qiw)が存在する。
これを再帰的に適用することで、a0x+b0y=dを満たす(x,y)の存在が言える。
また、はの倍数であるので、整数が存在してである。
すると、整数が存在しを満たすため、確かに整数解としてふさわしい。
そして整数解が1つあれば無数にあることを示したい。
これは難しくない。(x0,y0)が整数解であればax0+by0=cであるのでこれを変形してa(x0+db)+b(y0−da)=c。よって(x0+db,y0−da)が整数解になる。
これを好きなだけ繰り返すことでいくらでも異なる整数解を生成できる。よって整数解は無数にあると言える。
証明は省略するが、特殊解と整数を用いてと表せる解が、方程式の解のすべてである。
このことからはに整数解が存在するの中で絶対値最小であり、任意のはの倍数である、ということが言えます。
ということなので、をみたすをのベズー係数(Bézout coefficients)と名付けます(一意ではないことに気をつけてください)。
このベズー係数を1つ見つけることができるアルゴリズムとして知られているのが、拡張ユークリッドの互除法なのです!
最大公約数を求めるアルゴリズムであるユークリッドの互除法に機能を追加してベズー係数を求められるようにしていることから、この名前がつきました。
ここで驚くべき性質があります。拡張ユークリッドの互除法が与えるこのベズー係数はなんと、
をみたす2つのうち1つであるのです!
証明(クリックで開く)
帰納法を用いる。
一般性を失わずa>b>0とする。
まず拡張ユークリッドの互除法を実行してb=0に至る直前のステップ、つまりあるq≥1があってa=qd,b=dという場合を考える。
この場合与えられるベズー係数は(0,1)であるので、∣0∣≤∣∣∣dd∣∣∣,∣1∣≤∣∣∣∣dqd∣∣∣∣よりOK。
次にを満たすが条件を満たしていたと仮定した場合にを満たすが条件を満たすこと、すなわちとして
を示す。
∣x∣=∣w∣≤∣∣∣db∣∣∣より、∣x∣≤∣∣∣db∣∣∣
∣y∣=∣z−qw∣≤∣z∣+∣qw∣=∣∣∣dr∣∣∣+∣∣∣∣dqb∣∣∣∣=∣∣∣∣dqb+r∣∣∣∣(∵bq>0,r>0)=∣∣∣da∣∣∣
より、∣y∣≤∣∣∣da∣∣∣
以上より、機能法を用いて拡張ユークリッドの互除法が与えるベズー係数が
を満たすものであることが示された。
2項のベズー係数を拡張ユークリッドの互除法を用いて計算できること、わかっていただけましたか?
C++による実装例はこちらから見れます。
3項でのベズー係数
ここからが本題です
N項での議論をする前に前準備として3項での議論をしようと思います。
実はベズー係数は3項としてもその存在が言えます。つまり、を0でない整数として
をみたす整数(x,y,z)が存在します(そしてはそれを存在たらしめる右項のうち絶対値最小)。
これは次のようにして存在証明・計算ができます。
- のベズー係数を拡張ユークリッドの互除法を用いて一つ求める。
- のベズー係数を拡張ユークリッドの互除法を用いて一つ求める。
- がのベズー係数の一つである。
これが何故かというと、
の定義より、とが成り立ち、2つ目の式の最初の部分に1つ目の式を代入することで、
を導けるからです。
なので、はベズー係数の定義を満たしています。
N項でのベズー係数
書いてて思ったんだけどニッチすぎませんか?競プロですら出てくるか怪しいぞ
3項のものと同様の手順でN項に拡張できます。
最後の1項を除くN-1項のベズー係数が見つかっているものとします。
つまり、がのベズー係数の一つであることがわかっているものとします。また、とします。
- のベズー係数を拡張ユークリッドの互除法を用いて一つ求める
- がのベズー係数の一つである。
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
気休め程度が気に食わないので、オーダーレベルで対策をします。
を合体すると、
を得ます。
つまり、のベズー係数とのベズー係数がわかっていれば、
- のベズー係数を拡張ユークリッドの互除法を用いて一つ求める。
- がのベズー係数の一つである。
という手順でのベズー係数がで計算できます。
ではとをどうやって計算するかというと、再帰を使います。
セグ木や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と呼ばれるテクニックである、
互いに素な整数と整数があり、
であるとき
を満たす整数を求める。
というアルゴリズムは、拡張ユークリッドの互除法を用いて実行できます!
- のベズー係数を拡張ユークリッドの互除法を用いて一つ求める
- (をで割ったあまり)である。
あとがき
最初は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