元祖ワシ的日記

一つ前のエントリ

法Nでの逆数を求める

の手法は途中の計算結果を保存しておかないといけないという問題がありました。そこでそれを改良したアルゴリズムを紹介します。このアルゴリズムもユークリッドの互除法を使うという点では同じです。違うのは

\[ax + by = 1\]

ではなく

\[ax + by = N\]

を考えている点がです。この式は \[a=N, b=n\] とすると \[x=1, y=0\] という自明な解を持っています。

そこで \[a_1 = N, b_1 = n, x_1 = 1, y_1 = 0\] として前回と同じように

\[ a_{n+1} = b_n\\ b_{n+1} = m_n\\ x_{n+1} = d_n x_n + y_n\\ y_{n+1} = x_n \] と変形していきます。 \[ a_n x_n + b_n y_n = N \] が常に成り立っているので、

\[ a_n x_n = N - b_n y_n \] が成立します。 \[ a_n = b_{n-1}, y_n = x_{n-1} \] ですから \[ b_{n-1} x_n = N - b_n x_{n-1} \] と書き直すことができます。nを2からnまで変化させた等式

\[ b_1 x_2 = N - b_2 x_1\\ b_2 x_3 = N - b_3 x_2\\ \dots\\ b_{n-1} x_n = N - b_n x_{n-1} \]

について左辺どうし、右辺どうしを掛け合わせると

\[左辺 = b_1 b_2 \dots b_{n-1} x_2 x_3 \dots x_n \] \[右辺 = N(なんちゃら) \pm b_2 b_3 \dots b_{n} x_1 x_2 \dots x_{n-1} \]

となります。 $\pmod N$で考えれば $N(なんちゃら)$は0とみなせるので結局

\[b_1 x_n = \pm b_n x_1 \pmod N\]

が常に成立していることになります。符号はnの偶奇で決まります。なお $b_1, b_2, ...$ および $x_1, x_2, ...$がゼロでないという証明が別途必要ですがそれは最後に記します。

ところで、$b_n$はいずれ1になります。$x_1 = 1, b_1 = n$でしたから、その時点で

\[n x_n = \pm 1\]

となり $x_n$が求める解になります。このアルゴリズムは途中の計算結果を保存しておかなくてよいという利点があり一つ前のエントリより優れていると言えましょう。

サンプルコード

#!/usr/bin/env ruby
    
n = a = 916132831
b0 = b = 9973
x = 1
y = 0
sign = 1

while b != 1
  d = a / b
  m = a % b
  
  a = b
  b = m
  
  tmp = d * x + y
  y = x
  x = tmp
  
  sign = -sign
end

puts("#{sign * x} * #{b0} = #{b0 * x * sign % n} (mod #{n})")

さて、このアルゴリズムでも十分実用になるのですが、Nが奇数の場合はさらに効率よく求めることができます。 そのアルゴリズムを次のエントリで紹介しようと思います。

補足

$b_1, b_2, ...$ および $x_1, x_2, ...$がゼロでないという証明

$b_n$については単調に減少し1になったところでアルゴリズムは止まるので0でないことは明らか。 $x_n$については次の補題を利用します。

この補題より$\hbox{gcd}(x_n, y_n) = \hbox{gcd}(x_1, y_1) = 1$であり、したがって $x_n=0$ならば$y_n=1$でならねばならず それは$b_n=0$でなければならないので矛盾

と書くと何やら難しそうですが、簡単にいうならば$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では使われています。それを次回に紹介したいと思います。

三角関数使って表記すれば半径1の円に内接する正N角形の一辺の長さは

\[2\sin(\pi/N)\]

なんですが、この値は具体的に加減乗除と累乗根を使って表すことができます。 それを計算してみようという試み。

正三角形

\[2\sin(\pi/3) = \sqrt{3}\]

正方形

\[2\sin(\pi/4) = \sqrt{2}\]

正6角形

\[2\sin(\pi/6) = 1\]

ここまでは図を書けば殆ど自明です。正8角形、正12角形は半角の公式を使うと比較的簡単に求まります。半角の公式とは

\[\sin^2(\alpha/2) = \frac{1−\cos\alpha}{2}\]

です。

正8角形

\begin{eqnarray*} 2\sin(\pi/8) = 2\sqrt{\frac{1-\cos(\pi/4)}{2}}\\ = \sqrt{2-\sqrt{2}} \end{eqnarray*}

正12角形

\begin{eqnarray*} 2\sin(\pi/12) = 2\sqrt{\frac{1-\cos(\pi/6)}{2}}\\ = \sqrt{2-\sqrt{3}} \end{eqnarray*}

根号の中に根号になりますね。この先、16, 32, 64または24, 48, ..と増えていっても、半角の公式でどんどん数が求まっていきます。もちろん増える度に根号は深くなります、こんな風に辺の数が増えていくと、根号がどんどん深くなると予想されるのですが、正10角形は比較的シンプルな値になります。

正10角形

\[2\sin(\pi/10) = \frac{-1+\sqrt{5}}{2}\]

となります。正5角形よりシンプル。 実は正10角形は正5角形の作図がどうやって行われるかが理解できていると簡単に求めることができます。

一般に$x^n-1$の根を複素平面にプロットすると単位円に内接する正n角形の頂点となります。正10角形をプロットすると図のようになり、 一辺の長さは図の$P_1Q$になりますが、 これが$P_1+P_4$と原点Oからの距離($|P_1+P_4|$)に等しいことが見て取れると思います。

$P_1, P_2, P_3, P_4$は正5角形の頂点、すなわち$x^5-1$の根です。これには1という自明な根があるので因数分解ができて、

\[ x^5-1 = (x-1)(x^4+x^3+x^2+x+1) \]

この第2項

\begin{eqnarray} x^4+x^3+x^2+x+1 \hbox{　　(★)} \end{eqnarray}

の根ということになります。その根の1つを$\zeta$とすれば、すべての根は

\[ \zeta, \zeta^2, \zeta^3, \zeta^4 \]

となります。また$\zeta$は★の根ですから$\zeta+\zeta^2+\zeta^3+\zeta^4 = -1$になります。 ここでは$P_1=\zeta$とします。 今

\begin{eqnarray*} X = \zeta + \zeta^4\\ Y = \zeta^2 + \zeta^3 \end{eqnarray*}

とすれば

\begin{eqnarray*} X+Y = \zeta+\zeta^2+\zeta^3+\zeta^4 = -1\\ XY = \zeta+\zeta^2+\zeta^3+\zeta^4 = -1 \end{eqnarray*}

ですから($\zeta^5 = 1$に注意）X, Yは

\[ x^2 + x - 1 \]

の根。よって

\[ X = \frac{-1 + \sqrt{5}}{2} \]

となります。Xが$P_1+P_4$なのは言わずもがなです。これから次の有名な問題を証明できます。

問題）円周率が3.05より大きいことを証明せよ

\[ 2\pi > 10 \times \frac{-1 + \sqrt{5}}{2} = 6.180... > 6.1 \]

模範解答では、正8角形や正12角形を使うものが多いわけですが、こちらの方がより根号の数が少なく計算が簡単ですね。$\sqrt{5}=2.236...$というのを知っていればいいわけで。

なお余談ながら

\[1: |P_2+P_3|\] \[|P_1+P_4|:1\]

はいわゆる黄金比(1:1.618....)です。こんな計算のやりかたを知らなくても正10角形には黄金比が隠れていることを知っていれば3.05より大きいのは直ぐにわかるでしょう。いきなり正10角形の辺の長さは黄金比1.618:1だから・・・と書き始めたら正解になるかどうかはわかりませんが。

\[ \alpha_{0} = \left(-1 + \sqrt{17}\right)/2 \] \[ \alpha_{1} = \left(-1 - \sqrt{17}\right)/2 \] \[ \beta_{0} = \left(\alpha_{0} + \sqrt{\alpha_{0}^2 +4}\right)/2 \] \[ \beta_{1} = \left(\alpha_{1} + \sqrt{\alpha_{1}^2 +4}\right)/2 \] \[ \gamma_{0} = \left(\beta_{0} + \sqrt{\beta_{0}^2 - 4\beta_{1}}\right)/2 \] \[ \delta_{0} = \left(\gamma_{0} + \sqrt{\gamma_{0}^2 - 4}\right)/2 \] \[ \delta_{0}^{17} = 1 \]

\[\int_a^b f(x) dx=F(b)-F(a)\]

\[\mathbf{V}_1 \times \mathbf{V}_2 = \begin{vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \\ \frac{\partial X}{\partial u} & \frac{\partial Y}{\partial u} & 0 \\ \frac{\partial X}{\partial v} & \frac{\partial Y}{\partial v} & 0 \end{vmatrix} \]