元祖ワシ的日記

と書くと何やら難しそうですが、簡単にいうならば$n$という整数が与えられたとき

$$nx = 1 \pmod N$$

となる数$x$を見つけようということです。このような逆数はnとNが互いに素なら必ずあります。法Nにおける逆数は整数の剰余の計算においてとても重要であり暗号の分野で頻繁に使われています。

(多くの）プログラム言語的に言い換えると

(n * x) % N = 1

となる数ということになるでしょう。ここで 細かいことをいうと$n$が負数のときどうなるん？という問題があります。例えばperlやrubyだと $n$が負の場合は（例えば -1 % 7)は 6を返しますが、CやNode.jsは-1を帰すようです。 ちなみに、数学的にはどちらでも問題なしです。何故どちらでもいいかを語り始めると長いので細かいことは有限体の本でも読んでもらうとして、結論から言うと

$$-1 = 6 \pmod 7$$

と思っておけば問題ないです。7を法とした世界では-1も6も同じものを指します。表現の仕方が違うだけです。

さてこのアプローチには大きく2種類あるようです。 一つはユークリッドの互除法を使うもの。もう一つはopensslの内部て使われているものです。後者のアルゴリズムに名前があるのかどうかはしりません。

最初にユークリッドの互除法を使う手法を紹介します。

ユークリッドの互除法とは数$a$, $b$が与えられたときにその最大公約数を求める手法ですが、その解を求める課程を利用して

$$ax + by = 1$$

を解くことができるのです。実際ここで $a=N, b=n$として両辺をNで割ったあまりをみれば

$$by = 1 \pmod N$$

ですからyが求めるものであることがわかります。ちなみに上の式は$a, b$が互いに素なら（つまりGCD(a,b)=1ならば必ず解を持つことが知られています。

さてそのアルゴリズムですが、このアルゴリズムをを一般式で書くと煩雑になるので具体的に $a = (N = )916132831, b = (n = ) 9973$ として説明します。 まず$a$を$b$で割った余りと商を利用して以下のように変形します。

$$916132831x + 9973y\\ = (91861\times 9973 + 3078)x + 9973y\\ = 3078x + 9973(91861x + y)$$

ここで $$x' = 91861x + y\\ y' = x$$ とすれば $$9973x' + 3078y'$$ という式が得られます。商と余りを利用してより小さい係数の式を得ることができました。

このように書くと何やら複雑ですがプログラム的には

a[n+1] = b[n]
b[n+1] = a[n] % b[n]

ということです。このように$a_n, b_n$だけを計算していくといつかは$a_n$は$1$に、$b_n$は0になります。そのとき

$$x_n = 1\\ y_n = 0$$ とすると、 $$a_n x_n + b_n y_n = 1$$ を満たしますから、 これから逆に、$x_{n-1}, y_{n-1}$が求まってゆき、最終的に$y_1$が求まります。 これが求める答えになります。 例として、$a_1=916132831$と$b_1=9973$で計算していくと以下のようになります。

         #         a         b         d         m
         1 916132831      9973     91861      3078
         2      9973      3078         3       739
         3      3078       739         4       122
         4       739       122         6         7
         5       122         7        17         3
         6         7         3         2         1
         7         3         1         3         0
         8         1         0

これから$x_8=1, y_8=0$が求まります。 これから$x_7, y_7, x_6, y_6, ...$と計算していくと

         #         x         y
         7         0         1
         6         1        -2
         5        -2        35
         4        35      -212
         3      -212       883
         2       883     -2861
         1     -2861 262815204

となり、$y_1 = 262815204$が求まります。

実際 $9973 \times 262815204 = 1 \pmod N$です。

サンプルコード

#!/usr/bin/env ruby

a = []
b = []
d = []
m = []
x = []
y = []
n = a[1] = 916132831
b[1] = 9973

i = 1
while b[i] != 0
  d[i] = a[i] / b[i]
  m[i] = a[i] % b[i]

  a[i + 1] = b[i]
  b[i + 1] = m[i]
  
  i = i + 1
end

x[i] = 1
y[i] = 0
j = i
while j > 1
  j = j - 1
  x[j] = y[j + 1]
  y[j] = x[j + 1] - d[j] * x[j]
end

puts("#{y[1]} * #{b[1]} = #{b[1] * y[1] % n} (mod #{n})") 

ユークリッドの互除法はとても早く1まで係数が降下するので効率もよいのですが$a_n, b_n$の途中の計算結果を保存しておかないといけないという問題があります。また、$d, m$を求めるのに商計算をしなくてはなりません。 一般的に商の計算はとても遅く、特にNが数1000bitともなるとその計算がオーバヘッドとなるので、出来ることなら使いたくないというのが本音です。そこでより効率のよいアルゴリズムがopensslでは使われています。それを次回に紹介したいと思います。