内容は剰余の乗法の逆元(モジュラ逆数)についてです。
CTFをやると、RSA暗号の秘密鍵を解く問題で、次のようなWriteupを見たことがあるかと思います。
import libnum
e = 65537
p = 12476682960795779723419989287306239606331347310604553825605263028855086418051086300006278888049896375754096827163121306696417314531666670662341673511789487
q = 10626608485185909739187126602183513204215955466032868268570782358820047751795982373252873333276843831383542360202874683262678664975380865415068554992090483
c = 126028558760741438230925566962334702896791270808414391828894437291120207835997354509410869568173448979784329516392078311028442458946906210498176026685581176779209904039179175088984829699783861789249260322822491502790157863487861171337660785254139478229397872969136220001367017765809130698515063339265165085655
N = p* q
phi = ( p- 1 ) * ( q- 1 )
d = libnum. modular. invmod( e, phi)
print ( libnum. n2s( pow ( c, d, N) ) )
いわゆる、N N N が素因数分解で p × q p \times q p × q にできたら秘密鍵 d d d は計算できるというやつですが、ここで出てくるメソッド libnum.modular.invmod()
は数学では何をしているか説明したいと思います。
まず結論からいうと、invmod(e, phi)
は次のように書き換えることができます。
d ≡ i n v m o d ( e , ϕ ) d ≡ e − 1 ( m o d ϕ ) \begin{aligned}
d &\equiv \mathrm{invmod}(e, \phi) \\
d &\equiv e^{-1} \pmod{\phi}
\end{aligned} d d ≡ i n v m o d ( e , ϕ ) ≡ e − 1 ( m o d ϕ )
ただし、d d d , e e e , ϕ \phi ϕ は整数です。
ここで剰余演算について詳しく知らない人は「逆数の剰余って計算できるの?」と疑問に思います。
おそらく疑問に思っているところは次の点だと思います。
例えば e = 7 , ϕ = 5 e = 7, \phi = 5 e = 7 , ϕ = 5 として、
\gdef\mod{ {\;\mathrm{mod}\;} }
7 7 7 の逆数は 7 − 1 = 1 7 7^{-1} = \dfrac{1}{7} 7 − 1 = 7 1 (有理数)
7 7 7 を5 5 5 で割った余りは 7 m o d 5 7 \mod 5 7 m o d 5 で答えは 2 2 2 (剰余演算は整数のみ)
なので、7 − 1 ( m o d 5 ) 7^{-1} \pmod{5} 7 − 1 ( m o d 5 ) は一見すると有理数の剰余演算に見えるので、計算できないのではないかと思ってしまいます。
1. 剰余と合同式
まず、数学のおさらいをしていきたいと思います。
C言語などでは %
演算子で剰余を求めることができますが、改めて剰余の定義を数式で定義しておきましょう。
余りの定義
a a a を b b b で割ったときの商 q q q と余り r r r を次のように表す。ただし、a , b a,b a , b は自然数、q , r q,r q , r は自然数または 0 0 0 とする。
a = b q + r ( 0 ≤ r < b ) (1.1) a = bq + r \;\;\;(0 \le r < b)
\tag{1.1} a = b q + r ( 0 ≤ r < b ) ( 1 . 1 )
剰余の定義
整数 a , b , q , r a,b,q,r a , b , q , r で、b ≠ 0 b \ne 0 b = 0 とする。
a m o d b = r ⟺ a = b q + r ( 0 ≤ r < ∣ b ∣ ) (1.2) a \mod b = r \;\Longleftrightarrow\; a = bq + r \;\;\;(0 \le r < |b|)
\tag{1.2} a m o d b = r ⟺ a = b q + r ( 0 ≤ r < ∣ b ∣ ) ( 1 . 2 )
要するに、自然数の余りの定義を整数に広げたものが剰余(m o d \mathrm{mod} m o d )です。
例えば 7 m o d ( − 3 ) 7 \mod (-3) 7 m o d ( − 3 ) の答えは何になるでしょうか?
式(1.2)に沿って置き換えると以下のようになります。
7 m o d ( − 3 ) = r ⟺ 7 = ( − 3 ) q + r ( 0 ≤ r < ∣ − 3 ∣ ) 7 \mod (-3) = r \Longleftrightarrow 7 = (-3) q + r \;\;\;(0 \le r < \lvert -3\rvert) 7 m o d ( − 3 ) = r ⟺ 7 = ( − 3 ) q + r ( 0 ≤ r < ∣ − 3 ∣ )
条件 0 ≤ r < 3 0 \le r < 3 0 ≤ r < 3 の下で、整数 q = − 2 , r = 1 q = -2, r = 1 q = − 2 , r = 1 とすると 7 = ( − 3 ) ( − 2 ) + 1 7 = (-3)(-2) + 1 7 = ( − 3 ) ( − 2 ) + 1 となり上式を満たすので、答えは 7 m o d ( − 3 ) = 1 7 \mod (-3) = 1 7 m o d ( − 3 ) = 1 となります。
余談ですが、プログラミング言語によって剰余の結果は異なります。
剰余の次は合同式について説明します。
合同式の定義
a , b a,b a , b を c c c で割った余りが等しいとき「a , b a,b a , b は c c c を法 として合同 」といい、次のように表す。
a ≡ b ( m o d c ) (1.3) a \equiv b \pmod{c}
\tag{1.3} a ≡ b ( m o d c ) ( 1 . 3 )
合同式の加算・減算・乗算は一般的な等式と同じです。
a + C ≡ b + C ( m o d N ) a − C ≡ b − C ( m o d N ) a × C ≡ b × C ( m o d N ) \begin{aligned}
a + C &\equiv b + C \pmod N \\
a - C &\equiv b - C \pmod N \\
a \times C &\equiv b \times C \pmod N
\end{aligned} a + C a − C a × C ≡ b + C ( m o d N ) ≡ b − C ( m o d N ) ≡ b × C ( m o d N )
ただし、合同式の除算はできる場合とできない場合があります。
除算ができる条件は、割る数 C C C と法 N N N が互いに素であること (C C C と N N N の最大公約数が1 1 1 であること)です。
a ÷ C ≡ b ÷ C ( m o d N ) ( C ⊥ N ) a \div C \equiv b \div C \pmod N \;\;\;\;\; (C \bot N) a ÷ C ≡ b ÷ C ( m o d N ) ( C ⊥ N )
除算できない例(割る数 3 3 3 と法 6 6 6 が互いに素ではない):
9 ≡ 21 ( m o d 6 ) 3 ≢ 7 ( m o d 6 ) \begin{aligned}
9 &\equiv 21 \pmod 6 \\
3 &\not\equiv 7 \pmod 6
\end{aligned} 9 3 ≡ 2 1 ( m o d 6 ) ≡ 7 ( m o d 6 )
除算できる例(割る数 3 3 3 と法 5 5 5 が互いに素):
51 ≡ 21 ( m o d 5 ) 17 ≡ 7 ( m o d 5 ) \begin{aligned}
51 &\equiv 21 \pmod{5} \\
17 &\equiv 7 \pmod{5}
\end{aligned} 5 1 1 7 ≡ 2 1 ( m o d 5 ) ≡ 7 ( m o d 5 )
この「割る数」と「法」が互いに素であるという条件は、合同式の除算が可能になるだけではなくて、乗算の逆元の存在を示す十分条件にもなります。
2. 剰余の乗算の逆元
ここからは「割り算をする」を「逆数を掛け算する」と捉えることにします。
そこでまず、乗算の逆元を理解するには単位元と逆元について理解する必要があります。
なお、集合の「元」とは集合の「要素」と同じ意味です。
例えば、整数全体の集合 Z \mathbb{Z} Z の加算における単位元は 0 0 0 で、乗算における単位元は 1 1 1 です。
例えば、整数全体の集合 Z \mathbb{Z} Z の加算では3 3 3 の逆元は− 3 -3 − 3 になりますが、Z \mathbb{Z} Z の乗算では 4 × b = 1 4 \times b = 1 4 × b = 1 を満たす整数 b b b は存在しないので、乗算の逆元は存在しません。
集合を有理数全体の集合 Q \mathbb{Q} Q にすると4 4 4 の乗算の逆元は存在し、4 − 1 = 1 4 4^{-1} = \frac{1}{4} 4 − 1 = 4 1 となります。
ここで、剰余の( m o d N ) \pmod{N} ( m o d N ) の世界でも逆元が定義できれば一番最初の疑問を解くことができそうです。
つまり 7 × b ≡ 1 ( m o d 5 ) 7 \times b \equiv 1 \pmod 5 7 × b ≡ 1 ( m o d 5 ) を満たす整数 b b b があるとすれば、b ≡ 7 − 1 ( m o d 5 ) b \equiv 7^{-1} \pmod 5 b ≡ 7 − 1 ( m o d 5 ) となるので、b b b は 7 7 7 の逆元であると言えます。
なお、( m o d N ) \pmod{N} ( m o d N ) の世界を集合で表すときは、この集合を剰余環 と呼び、
Z / N Z = { 0 , 1 , 2 , . . . , N − 1 } \mathbb{Z}/N\mathbb{Z} = \{0,1,2,...,N-1\} Z / N Z = { 0 , 1 , 2 , . . . , N − 1 }
と書きます。例えば、( m o d 5 ) \pmod{5} ( m o d 5 ) の世界は集合で次のように表せます。
Z / 5 Z = { 0 , 1 , 2 , 3 , 4 } \mathbb{Z}/5\mathbb{Z} = \{0,1,2,3,4\} Z / 5 Z = { 0 , 1 , 2 , 3 , 4 }
剰余環は環の性質を持っています。なので、剰余の加算・乗算は整数の加算・乗算とほとんど変わりません。
では、集合 Z / 5 Z \mathbb{Z}/5\mathbb{Z} Z / 5 Z における 2 ( m o d 5 ) 2\pmod{5} 2 ( m o d 5 ) の逆元である 2 − 1 ( m o d 5 ) 2^{-1}\pmod{5} 2 − 1 ( m o d 5 ) とは一体何でしょうか?
集合は Z / 5 Z = { 0 , 1 , 2 , 3 , 4 } \mathbb{Z}/5\mathbb{Z} = \{0,1,2,3,4\} Z / 5 Z = { 0 , 1 , 2 , 3 , 4 } なので、全部の値をそれぞれ掛け合わせてみて乗法の単位元の 1 1 1 になるか計算してみましょう。
( 2 ( m o d 5 ) ) × ( 0 ( m o d 5 ) ) = 0 ( m o d 5 ) ≡ 0 ( 2 ( m o d 5 ) ) × ( 1 ( m o d 5 ) ) = 2 ( m o d 5 ) ≡ 2 ( 2 ( m o d 5 ) ) × ( 2 ( m o d 5 ) ) = 4 ( m o d 5 ) ≡ 4 ( 2 ( m o d 5 ) ) × ( 3 ( m o d 5 ) ) = 6 ( m o d 5 ) ≡ 1 ( 2 ( m o d 5 ) ) × ( 4 ( m o d 5 ) ) = 8 ( m o d 5 ) ≡ 3 \begin{aligned}
(2 \;(\mathrm{mod}\; 5)) \times (0 \;(\mathrm{mod}\; 5)) \;\;\;
&= 0 \;(\mathrm{mod}\; 5)
&\equiv 0 \\
(2 \;(\mathrm{mod}\; 5)) \times (1 \;(\mathrm{mod}\; 5)) \;\;\;
&= 2 \;(\mathrm{mod}\; 5)
&\equiv 2 \\
(2 \;(\mathrm{mod}\; 5)) \times (2 \;(\mathrm{mod}\; 5)) \;\;\;
&= 4 \;(\mathrm{mod}\; 5)
&\equiv 4 \\
(2 \;(\mathrm{mod}\; 5)) \times (3 \;(\mathrm{mod}\; 5)) \;\;\;
&= 6 \;(\mathrm{mod}\; 5)
&\equiv 1 \\
(2 \;(\mathrm{mod}\; 5)) \times (4 \;(\mathrm{mod}\; 5)) \;\;\;
&= 8 \;(\mathrm{mod}\; 5)
&\equiv 3 \\
\end{aligned} ( 2 ( m o d 5 ) ) × ( 0 ( m o d 5 ) ) ( 2 ( m o d 5 ) ) × ( 1 ( m o d 5 ) ) ( 2 ( m o d 5 ) ) × ( 2 ( m o d 5 ) ) ( 2 ( m o d 5 ) ) × ( 3 ( m o d 5 ) ) ( 2 ( m o d 5 ) ) × ( 4 ( m o d 5 ) ) = 0 ( m o d 5 ) = 2 ( m o d 5 ) = 4 ( m o d 5 ) = 6 ( m o d 5 ) = 8 ( m o d 5 ) ≡ 0 ≡ 2 ≡ 4 ≡ 1 ≡ 3
Z / 5 Z \mathbb{Z}/5\mathbb{Z} Z / 5 Z では2 2 2 に3 3 3 を掛けると乗法の単位元 1 1 1 となることから、
2 ( m o d 5 ) 2 \pmod 5 2 ( m o d 5 ) の逆元は 3 ( m o d 5 ) 3 \pmod 5 3 ( m o d 5 ) であることがわかりました。
2 − 1 ( m o d 5 ) = 3 ( m o d 5 ) ≡ 3 \begin{aligned}
2^{-1} \;(\mathrm{mod}\; 5) &= 3 \;(\mathrm{mod}\; 5) \\
&\equiv 3
\end{aligned} 2 − 1 ( m o d 5 ) = 3 ( m o d 5 ) ≡ 3
2 2 2 と7 7 7 は5 5 5 を法として合同なので、2 2 2 と7 7 7 の逆元も同じになります。
7 − 1 ( m o d 5 ) = ( 7 ( m o d 5 ) ) − 1 = ( 2 ( m o d 5 ) ) − 1 = 2 − 1 ( m o d 5 ) = 3 ( m o d 5 ) ≡ 3 \begin{aligned}
7^{-1} \;(\mathrm{mod}\; 5)
&= (7 \;(\mathrm{mod}\; 5))^{-1} \\
&= (2 \;(\mathrm{mod}\; 5))^{-1} \\
&= 2^{-1} \;(\mathrm{mod}\; 5) \\
&= 3 \;(\mathrm{mod}\; 5) \\
&\equiv 3
\end{aligned} 7 − 1 ( m o d 5 ) = ( 7 ( m o d 5 ) ) − 1 = ( 2 ( m o d 5 ) ) − 1 = 2 − 1 ( m o d 5 ) = 3 ( m o d 5 ) ≡ 3
では、0 0 0 は何を掛けても0 0 0 になるので除くとして、集合 Z / 5 Z \mathbb{Z}/5\mathbb{Z} Z / 5 Z から0 0 0 を除いた集合 ( Z / 5 Z ) ∗ (\mathbb{Z}/5\mathbb{Z})^* ( Z / 5 Z ) ∗ の任意の元について、乗算の逆元は存在すると言えるのでしょうか?
Z / 4 Z \mathbb{Z}/4\mathbb{Z} Z / 4 Z と Z / 5 Z \mathbb{Z}/5\mathbb{Z} Z / 5 Z についてそれぞれの乗算の結果をまとめた表を以下に示します。この表のことを演算表 と呼び、乗算についての演算表を乗算表と呼んだりします。
逆元が存在するかを調べるために、乗算の結果が1なら赤色にしている
そうすると、乗算表からわかるように ( Z / 5 Z ) ∗ (\mathbb{Z}/5\mathbb{Z})^* ( Z / 5 Z ) ∗ の全ての元には逆元が存在する(つまり乗算したら1 1 1 になる元が存在する)ことが確認できますが、( Z / 4 Z ) ∗ (\mathbb{Z}/4\mathbb{Z})^* ( Z / 4 Z ) ∗ の元である 2 2 2 は任意の元と掛け合わせても1 1 1 にならないので逆元は存在しないことになります。
それでは、どのような条件を満たす時に ( Z / N Z ) ∗ (\mathbb{Z}/N\mathbb{Z})^* ( Z / N Z ) ∗ の全ての元に逆元が存在すると言えるのでしょうか?
ここで合同式の除算ができる条件を思い出すと「割る数 C C C と法 N N N が互いに素であること」でした。
つまり、任意の元で除算ができる、すなわち逆元が存在することを示すには、
集合 ( Z / N Z ) ∗ = { 1 , 2 , … , N − 1 } (\mathbb{Z}/N\mathbb{Z})^* = \{1,2,\ldots,N-1\} ( Z / N Z ) ∗ = { 1 , 2 , … , N − 1 } のそれぞれの元 1 , 2 , … , N − 1 1,2,\ldots,N-1 1 , 2 , … , N − 1 が法 N N N と互いに素であることが条件になります。
この条件を満たす法 N N N は素数しかないのです。
p p p を素数とすると、剰余環 Z / p Z \mathbb{Z}/p\mathbb{Z} Z / p Z の全ての元について逆元が存在します。
この集合を体 といい、素数 p p p を法とする剰余環 Z / p Z \mathbb{Z}/p\mathbb{Z} Z / p Z を有限体 と呼びます。
F p = Z / p Z ( p is a prime ) \mathbb{F}_p = \mathbb{Z}/p\mathbb{Z} \;\;\;\;\; (p \;\text{is a prime}) F p = Z / p Z ( p is a prime )
ちなみに、RSA暗号に関して言えば有限体は全く出てこないので知らなくても問題ないですが、せっかく剰余環について書いたので有限体にも触れておこうかなと思った次第です。
ここまでで剰余の乗算の逆元について説明してきたので、次はいよいよRSAの話をしたいと思います。
3. RSA暗号方式の鍵生成
逆元を求めるために毎回乗算表を作っていたら大変ですよね。
実際、秘密鍵を作るときは大きな数(例えば600桁)を使うので、乗算表を作っていたら日が暮れてしまいます。
では、RSAの鍵生成ではどのように秘密鍵を作っているか確認してみましょう。
RSA暗号のアルゴリズムは次のように定義されています。
RSA:鍵生成
二つの大きな素数 p , q p, q p , q を適切に決める。
N = p q N = pq N = p q を求める。
オイラーの関数 φ ( N ) = ( p − 1 ) ( q − 1 ) \varphi(N) = (p-1)(q-1) φ ( N ) = ( p − 1 ) ( q − 1 ) を計算する。
整数 e e e を φ ( N ) \varphi(N) φ ( N ) 未満の正の整数で、 φ ( N ) \varphi(N) φ ( N ) と互いに素な数とする。
公開鍵を n , e n, e n , e とする。
秘密鍵 d d d を、φ ( N ) \varphi(N) φ ( N ) を法とした e e e の逆数(e d ≡ 1 ( m o d φ ( N ) ) ed \equiv 1 \pmod{\varphi(N)} e d ≡ 1 ( m o d φ ( N ) ) )とする。d d d は拡張ユークリッドの互除法を使えば容易に求まる。
RSA:暗号化
m m m を平文空間 Z n = { 0 , 1 , 2 , … , n − 1 } \mathbb{Z}_n = \{0,1,2,\ldots,n-1\} Z n = { 0 , 1 , 2 , … , n − 1 } の元とすると、暗号文 c c c は次の式で求まる。
c = m e m o d n (3.1) c = m^e \mod n
\tag{3.1} c = m e m o d n ( 3 . 1 )
RSA:復号
c c c を暗号文とすると、平文 m m m は次の式で求まる。
m = c d m o d n (3.2) m = c^d \mod n
\tag{3.2} m = c d m o d n ( 3 . 2 )
鍵生成のところで e d ≡ 1 ( m o d φ ( N ) ) ed \equiv 1 \pmod{\varphi(N)} e d ≡ 1 ( m o d φ ( N ) ) と書かれているのは e ≡ d − 1 ( m o d φ ( N ) ) e \equiv d^{-1} \pmod{\varphi(N)} e ≡ d − 1 ( m o d φ ( N ) ) と同じ式です。
φ ( N ) \varphi(N) φ ( N ) とはオイラーの関数で、1 1 1 から N N N までの自然数のうち N N N と互いに素なものの個数を返す関数のことです。例えば、N N N が素数なら 1 1 1 から N − 1 N-1 N − 1 までの自然数の全てと互いに素なので、φ ( N ) = N − 1 \varphi(N) = N-1 φ ( N ) = N − 1 となります。また、オイラーの関数は互いに素な数の積に関して乗法的であるので、N N N が二つの素数の合成数 N = p q N = pq N = p q の場合は φ ( N ) = ( p − 1 ) ( q − 1 ) \varphi(N) = (p-1)(q-1) φ ( N ) = ( p − 1 ) ( q − 1 ) が成り立ちます。
RSA暗号の安全性は素因数分解の困難性に依拠しています。「依拠する」とは暗号の解読それらの問題の困難さとの間に厳密な等価性が成り立つという意味ではなく、ただ頼っているという意味です。
N N N の素因数を知らない解読者は、法 φ ( N ) \varphi(N) φ ( N ) を知らないので、ユークリッドの互除法を使えず、秘密鍵は求められません。しかし今後コンピュータの性能が向上して N N N の素因数分解が現実的な時間でできるようになったとき、簡単に N N N は素因数分解されて秘密鍵 d d d は計算で求められてしまっては困ります。
なので、今使っているRSAのビット数が心配という方は、CRYPTRECが公表している『CRYPTREC暗号リスト(電子政府推奨暗号リスト)』に一度目を通しておくと良いです。
RSAの暗号化して復号したときの結果が等しくなることはオイラーの定理 a φ ( n ) ≡ 1 ( m o d n ) a^{\varphi(n)} \equiv 1 \pmod{n} a φ ( n ) ≡ 1 ( m o d n ) を使うと簡単に説明できます(ただし、m < n m < n m < n とする)。
m = c d m o d n = ( m e m o d n ) d m o d n = m e d m o d n = m 1 + k φ ( n ) m o d n = m 1 ⋅ m k φ ( n ) m o d n = m ⋅ 1 m o d n = m \begin{aligned}
m &= c^d \mod n = (m^e \mod n)^d \mod n \\
&= m^{ed} \mod n = m^{1 + k \varphi(n)} \mod n \\
&= m^1 \cdot{} m^{k\varphi(n)} \mod n = m \cdot{} 1 \mod n = m
\end{aligned} m = c d m o d n = ( m e m o d n ) d m o d n = m e d m o d n = m 1 + k φ ( n ) m o d n = m 1 ⋅ m k φ ( n ) m o d n = m ⋅ 1 m o d n = m
話を戻して、秘密鍵 d d d は拡張ユークリッドの互除法を使うと(乗算表を使わなくても)簡単に求めることができます。それでは次は拡張ユークリッドの互除法の話でもしたいと思います。
4. 拡張ユークリッドの互除法による invmod
冒頭に示したメソッド invmod
は拡張ユークリッドの互除法を使うことで簡単に求めることができます。
そこで、ユークリッドの互除法と拡張ユークリッドの互除法について説明したいと思います。
4.1. ユークリッドの互除法
拡張ユークリッドの互除法の前に、ユークリッドの互除法の説明を簡単にしておきます。
ユークリッドの互除法とは最大公約数を求めるアルゴリズムです。
余りの定義(1.1)より、a a a を b b b で割った商 q q q と余り r r r とすると、a a a は a = b q + r a = bq + r a = b q + r のように表すことができます。
r = 0 r = 0 r = 0 のときは明らかに a a a と b b b の最大公約数 gcd ( a , b ) \gcd(a,b) g cd( a , b ) は b b b となります。
r ≠ 0 r \ne 0 r = 0 のときは r = a − b q r = a - bq r = a − b q と表せるので、gcd ( a , b ) \gcd(a,b) g cd( a , b ) は r r r の約数です。なので gcd ( b , r ) ≤ gcd ( a , b ) \gcd(b,r) \le \gcd(a,b) g cd( b , r ) ≤ g cd( a , b ) となります。
一方、a = b q + r a = bq + r a = b q + r と表されるから、a a a は gcd ( b , r ) \gcd(b,r) g cd( b , r ) の約数です。なので gcd ( b , r ) ≥ gcd ( a , b ) \gcd(b,r) \ge \gcd(a, b) g cd( b , r ) ≥ g cd( a , b ) となります。
したがって、gcd ( b , r ) = gcd ( a , b ) \gcd(b,r) = \gcd(a,b) g cd( b , r ) = g cd( a , b ) となるので、gcd ( b , r ) \gcd(b,r) g cd( b , r ) を求めることは gcd ( a , b ) \gcd(a,b) g cd( a , b ) を求めることと同じということになります。
( a , b ) (a,b) ( a , b ) よりも ( b , r ) (b,r) ( b , r ) の方が小さい整数の組なので、当然 gcd ( a , b ) \gcd(a,b) g cd( a , b ) よりも gcd ( b , r ) \gcd(b,r) g cd( b , r ) を求める方が簡単です。説明を端折りますが、この処理を繰り返していって r = 0 r = 0 r = 0 となったときの b b b が求めるべき最大公約数となります。
最大公約数を求めるプログラムをPythonで書くと下のようになります。
def gcd ( a, b) :
while b != 0 :
a, b = b, a % b
return a
4.2. 拡張ユークリッドの互除法
ユークリッドの互除法の応用として、拡張ユークリッドの互除法というアルゴリズムがあります。
拡張ユークリッドの互除法は a , b a,b a , b を正の整数とすると、次のベズーの等式を満たす適当な整数 x , y x,y x , y が存在するので、この x , y x,y x , y を求めるアルゴリズムです。
a x + b y = gcd ( a , b ) ax + by = \gcd(a,b) a x + b y = g cd( a , b )
ユークリッドの互除法において、その k k k 回目の除算の商を q k q_k q k 、余りを r k r_k r k とすると、ユークリッドの互除法は以下のように書くことができます。ただし r k = 0 r_k = 0 r k = 0 になったときにユークリッドの互除法は終了します。
r 1 = a − q 1 b r 2 = b − q 2 r 1 r 3 = r 1 − q 3 r 2 r 4 = r 2 − q 4 r 3 ⋮ r k = r k − 1 − q k r k − 1 ⋮ \begin{aligned}
r_1 &= a - q_1 b \\
r_2 &= b - q_2 r_1 \\
r_3 &= r_1 - q_3 r_2 \\
r_4 &= r_2 - q_4 r_3 \\
&\;\;\vdots \\
r_k &= r_{k-1} - q_k r_{k-1} \\
&\;\;\vdots
\end{aligned} r 1 r 2 r 3 r 4 r k = a − q 1 b = b − q 2 r 1 = r 1 − q 3 r 2 = r 2 − q 4 r 3 ⋮ = r k − 1 − q k r k − 1 ⋮
この式の右辺にある r k − 2 , r k − 1 r_{k-2}, r_{k-1} r k − 2 , r k − 1 を上の式で消去していくと、r k r_k r k は整数を係数とする a a a と b b b の線型結合で表すことができます。
例えば r 1 r_1 r 1 と r 2 r_2 r 2 について a a a と b b b の式に書き換えると次のようになります。
r 1 = 1 ⋅ a + ( − q 1 ) b r 2 = b − q 2 r 1 = ( − q 2 ) a + ( 1 + q 2 q 1 ) b \begin{aligned}
r_1 &= 1 \cdot{} a + (- q_1) b \\
r_2 &= b - q_2 r_1 = (-q_2)a + (1 + q_2 q_1)b \\
\end{aligned} r 1 r 2 = 1 ⋅ a + ( − q 1 ) b = b − q 2 r 1 = ( − q 2 ) a + ( 1 + q 2 q 1 ) b
r i = x i a + y i b r_i = x_i a + y_i b r i = x i a + y i b と表すことにすると、上式より x 1 = 1 , y 1 = − q 1 , x 2 = − q 2 , y 2 = 1 + q 2 q 1 x_1 = 1,\; y_1 = -q_1,\; x_2 = -q_2,\; y2 = 1 + q_2 q_1 x 1 = 1 , y 1 = − q 1 , x 2 = − q 2 , y 2 = 1 + q 2 q 1 となります。そこで、i < k i < k i < k を満たす i i i について r i r_i r i を表すと
r k = r k − 2 − q k r k − 1 = ( x k − 2 a + y k − 2 b ) − q k ( x k − 1 a + y k − 1 b ) = ( x k − 2 − q k x k − 1 ) a + ( y k − 2 − q k y k − 1 ) b = x k a + y k b \begin{aligned}
r_k &= r_{k-2} - q_k r_{k-1} \\
&= (x_{k-2} a + y_{k-2} b) - q_k (x_{k-1} a + y_{k-1} b) \\
&= (x_{k-2} - q_k x_{k-1}) a + (y_{k-2} - q_k y_{k-1}) b \\
&= x_k a + y_k b
\end{aligned} r k = r k − 2 − q k r k − 1 = ( x k − 2 a + y k − 2 b ) − q k ( x k − 1 a + y k − 1 b ) = ( x k − 2 − q k x k − 1 ) a + ( y k − 2 − q k y k − 1 ) b = x k a + y k b
より、x k , y k x_k, y_k x k , y k は次式になります。
x k = x k − 2 − q k x k − 1 y k = y k − 2 − q k y k − 1 \begin{aligned}
x_k &= x_{k-2} - q_k x_{k-1} \\
y_k &= y_{k-2} - q_k y_{k-1}
\end{aligned} x k y k = x k − 2 − q k x k − 1 = y k − 2 − q k y k − 1
初期条件を x 0 = 1 , y 0 = 0 , x 0 = 0 , x 1 = 1 x_0 = 1, y_0 = 0, x_0 = 0, x_1 = 1 x 0 = 1 , y 0 = 0 , x 0 = 0 , x 1 = 1 としてユークリッドの互除法でその都度求まる q k q_k q k を使って x i , y i x_i, y_i x i , y i を計算していけば、gcd ( a , b ) = r n = x n a + y n b \gcd(a,b) = r_n = x_n a + y_n b g cd( a , b ) = r n = x n a + y n b を満たす x n , y n x_n, y_n x n , y n を求めることができます。そして、このときの x n , y n x_n, y_n x n , y n が求めるべき x , y x, y x , y となります。
拡張ユークリッド互助法のアルゴリズムをPythonで書くと次の通りです。
この関数の返り値は (最大公約数gcd ( a , b ) \gcd(a,b) g cd( a , b ) の値, 整数x x x , 整数y y y ) です。
def xgcd ( a, b) :
x0, y0, x1, y1 = 1 , 0 , 0 , 1
while b != 0 :
q, a, b = a // b, b, a % b
x0, x1 = x1, x0 - q * x1
y0, y1 = y1, y0 - q * y1
return a, x0, y0
4.3. 剰余環の乗法逆元
剰余の乗法の逆元(invmod)は拡張ユークリッドの互除法を使って求めます。
もし a a a には法 m m m の乗法逆元が存在するなら、二つの整数 a , m a, m a , m の最大公約数 gcd ( a , m ) \gcd(a,m) g cd( a , m ) は 1 1 1 になります。すなわち、拡張ユークリッドの互除法で出てきたベズーの等式を使うと、次式が成り立ちます。
a x + m y = gcd ( a , m ) = 1 ax + my = \gcd(a,m) = 1 a x + m y = g cd( a , m ) = 1
書き換えると
a x = 1 + ( − y ) m ax = 1 + (-y)m a x = 1 + ( − y ) m
つまり
a x ≡ 1 ( m o d m ) ax \equiv 1 \pmod{m} a x ≡ 1 ( m o d m )
となるので、拡張ユークリッドの互除法を計算することで a a a の乗法逆元 x x x を求めることができます。
Pythonで書くと下のようになります。
def invmod ( a, m) :
g, x, y = xgcd( a, m)
if g != 1 :
raise Exception( 'modular inverse does not exist' )
else :
return x % m
試しに計算させてみましょう。
手で計算した 7 − 1 ( m o d 5 ) 7^{-1} \;(\mathrm{mod}\; 5) 7 − 1 ( m o d 5 ) や、Z / 4 Z \mathbb{Z}/4\mathbb{Z} Z / 4 Z の乗算表から逆元が存在しないことがわかっている 2 − 1 ( m o d 4 ) 2^{-1} \;(\mathrm{mod}\; 4) 2 − 1 ( m o d 4 ) を試してみます。
>>> invmod(7, 5)
3
>>> invmod(2, 4)
Exception: modular inverse does not exist
正しく計算できているようです。
という訳で invmod
もとい、剰余の乗法の逆元を求める話はこの辺で終わりにしたいと思います。
最後までお付き合い頂きありがとうございました。