-*-mode:text-*- * Publisher: 技術評論社 * Publish: SoftwareDesign * ThemeTitle: プログラムで理解する暗号 (cipher programming) * Auther: すずきひろのぶ ( Hirononou SUZUKI ) * E-Mail: hironobu -at- h2np -dot- net * Title: ソフトウェア・コンサルタント * TOPIC: RUBYのための道具 * 道具を揃える Rubyでの整数値計算では桁溢れを意識する必要がありません。整数の値はすべ てInegerクラスとして計算できます。32ビット以下の値の場合Fixnumというク ラスを使い、それ以上であれば多倍長整数のクラスBignumを自動的に使います。 もちろんFixnumクラスで計算していて、それで収まらなくなると自動的に Bignumクラスになります。a*x + b*y = ??? という計算のために一々桁数あふ れするかしないかプログラム上で気にすることもないわけです。 中身の話になりますがRubyの多倍長整数実装は、GNUの多倍長整数ライブラリ GMPといったものを使わず、多倍長整数を扱うアルゴリズムを使い独自にC言語 で書かれています。これはRubyの移植性向上のために、外部にあるオプショナ ルなライブラリを使わなかったのでしょう。そのため用意されているのは基本 的な演算のみです。 多倍長整数ライブラリGMPは、極限までスピードを追い求めるため、下位の演 算ライブラリはアセンブラで用意されています。AMD、Intel、SPARC、Alpha、 Motrolaは云うに及ばず、SH、VAXあるいはCRAYまで幅広くサポートしています。 また用意している関数も多く、ある数aの指数bの剰余mとした時のa^b mod m (modulo mでのa^b)を求めるmpz_powm()や、ある数a、剰余mとしたときの a^-1 mod m (modulo m でのaの逆数)を求めるmpz_invert()などがあり、数論 (Numeber Theory)の演算を行なう道具が揃っています。 GMPをRubyにリンクするのも1つの手ですが、それではつまらないので、Rubyの 道具を揃えて計算に備えましょう。ここでは練習も含めて最大公約数を求める 関数gcd(a,b)、modulo m でのaの逆数を求める関数invert(a,m)、そして高速 に大きな素数を判定する関数prime?(n)を作ります。 * 最大公約数gcd 最大公約数GCD(Greatest Common Divisor)を求める方法はこれはユークリッド の互除法という、基本中の基本であるアルゴリズムを使って解きます。aとbの 最大公約数を求める方法は次の通りです。 Step 0: a >0, b > 0, a >= bである aとbが与えられた時 Step 1: b != 0であるあいだ r = a mod b, a = b, b = rの計算を繰り返す。 Step 2: aが最大公約数GCDである。 ---Ruby--- def gcd(a,b) if ( a < b ); t=b; b=a; a=t; end while b > 0 r = a % b a = b ; b = r end return a end --------- なぜ、どこにでも用意されている関数gcdをわざわざ作ったかというと、これ は次のアルゴリズムを解説するためです。(筆者はアルゴリズムを思い出すた めです。) * modulo m でのa逆数 拡張ユークリッドの互除法について、簡単に説明しましょう。これはaとb、 d=gcd(a,b)である時 a*x + b*y = d を満たすxとyを求める時に利用する方法 です。この拡張ユークリッドの互除法に関してのアルゴリズムの説明は、Web 上にたくさんあります。 多くはDonald KnuthのThe Art of Computer Programming Vol2. pp.342 -- 343にある説明を基本としているようです。 d=gcd(a,b)である時 a*x + b*y = d を満たすxとyを求める Step 1: x1 = 0, x2 = 1, y1 = 1, y2 = 0とする Step 2: b > 0のあいだStep 3とStep 4を行なう Step 3: 次の計算を行なう q = a/b, r = a - p*b, x = x2 - p*x1, y = y2 - p*y1 Step 4: 次の代入を行なう a = b, b = r, x2 = x1, x1 = x, y2 = y1, y1 =y Step 5: d = a, x = x2, y = y2とする d, x, yが求める値である このアルゴリズムを使うのですが、これを応用して(簡略化して) gcd(a,m)=1 の時、a*x == 1 mod m 満たすxを求める、という問題に変形させます。 a*x+b*y=dの式でd=1、b=m、y=1であるような場合です。それをちょっと工夫を こらして短いプログラムにしました。プログラムを短くするために、けっこう 時間を食ったのでした...(こーゆーのにハマっちゃうよね)。 ところで、googleでは、拡張ユークリッドの互除法を解説したWebページはた くさん見つかりますが、modulo mでのaの逆数を計算するアルゴリズムのみを 書いたページが見つからなくて不思議な思いをしました。筆者が見つけたのは、 唯一、JAISTの宮地先生が筑波大学での行なわれた暗号の集中講義に使った資 料のなかの有限体上の除算についての章の応用問題として用意しているだけで した。暗号では頻繁に使うので、もうちょっと初心者解説や簡単なヒントがあっ てもいいような気がするのですが... 今、一生懸命、本棚を探しても見つからなかったのですが、このアルゴリズム と同じ方法の懇切親切な説明は、講談社ブルーバックの「情報セキュリティの 科学(太田・黒澤・渡辺)」にあったはずです。なぜかWebの上にはこのような 親切な方法が書いていません。 さて、プログラムですが最後にreturn (v+n) % n とせずに、return v % n と しているのは、Rubyの場合 vがマイナスの場合、(v+n) % nの値としてくれる ためです。そのためRuby以外の演算系ではv+nとした方が安全だと思います。 ---Ruby--- (変数名が上のアルゴリズム解説に対応していませんが、気にしないでください) def invert(x,n) a = n b = x % n p = 0 ; q = 1; v = 0; u = 1; while q > 0 p = a / b ; q = a % b w = v - (p*u) a = b ; b = q ; v = u ; u = w end return v % n end printf "22371 * (%d) mod 10093 = 1\n",invert(22371,10093) ---- Prof. Manindra and two of his students, Nitin Saxena and Neeraj * 大きな素数を高速に判定する 大きな数の素数判定には時間がかかる問題があります。中学校あたりでエラト ステネス(Eratosthenes)の篩(ふるい)、を習うと思いますし、コンピュータの 実習などでもこのアルゴリズムを実装するような練習問題が出されているので 一度ぐらいは書いたことがあるかも知れません。しかし、このようなアルゴリ ズム方法では、大きな素数には太刀打ちできません。あまりにも時間かかかり 過ぎるのです。最近、多項式時間での素数判定するAgarwal-Kayal-Saxenaアル ゴリズムが発表されましたがそれでもまだ実用的ではありません。 http://www.cse.iitk.ac.in/news/primality.html そのため素数かどうかを確率的に判定する高速なアルゴリズムを使います。こ のような判定方法をMonte Carlo法、かっこよく呼ぶ時はProbabilistic Prime Testingと呼びます。このような判定方法のなかで広く使われているのが Miller-Rabin法です。 ** Millerの判定法 「Nを整数とし、1 < a < N, N-1= 2^2 (s>=0,dは奇数)とする。a^(N-1) != 1 mod Nか、もしくは t = 2^r * dとし 1 < gcd(a^t -1,N) < Nであるr (0<=r 0 b[i] = (0x1 & v) v = v/2 i = i + 1 end return b end def miller_rabin(n,s) b=bitarray(n-1) i=b.size j =1 while j <= s a = 1 + (rand(n-2).to_i) if witness(a,n,i,b) == true return false end j+=1 end return true end def witness(a,n,i,b) d=1 x=0 while i > 0 x = d d = (d**2) % n if ( (d == 1) && (x != 1) && (x != (n-1)) ) return true end if ( b[i-1] == 1 ) d = (d * a ) % n end i -= 1 end if ( d != 1) return true end return false end ----- require "mathh" a= 2**1 + 1 while a < 100000 if ( a.prime? == true ) print a,"\n" end a += 2 end ----- * その他の道具 あとランダムな素数を生成する道具も用意しました。これはシステムが用意し ている疑似乱数デバイス/dev/urandomから必要なサイズのデータを読み込み、 それが素数かどうか判定しています。任意の疑似乱数デバイスに切替えること ができるようになっているので、/dev/urandomを信用せず、/dev/randomを使 うことも可能になっています。 整数においての割算での切りあげcdiv (ceil divide)は簡単にできます。Ruby のFloatクラスにはceilがあります。そこでceil{a/b}としたい時、 ((a*1.0)/b).ceilとaに1.0を乗算することでFloatに変換し、そこでceilを使っ て切り上げてしまうことができます。戻ってくるのは整数クラスです。 * おわりに これで公開鍵暗号などの計算に必要な基本的な道具は整いました。これは暗号 だけではなく数論で使う基本的な計算道具でもあります。なぜなら公開鍵暗号 は、数論の応用分野だからです。次回は実際に公開鍵暗号アルゴリズムや署名 アルゴリズムを実装してみましょう。 今回つくった道具は、下記URLに用意しておきます。まだ興味のある人は、使っ てみてください。 ----------------------------------------------------------------------- # # $Id: mathh.rb,v 1.5 2002/11/24 16:15:56 hironobu Exp $ # Simple Secure Stream Cryptography # Mathematical Library # Hironobu SUZUKI # LGPL, (C) 2002, Hironobu SUZUKI # Free free to send bug report to me. # class Random def initialize() @random_device="/dev/urandom" @verbose=false end def random_device=(str) @random_device=str end def verbose @verbose=true end def get(size) if (@verbose) ;STDERR.print "+";end value=0 r=File.open(@random_device) r.read(size).each_byte {|c| value = (value << 8) + c.to_i } r.close return value end def get_prime(lower,higher) bytesize = higher.bytesize while true r = self.get(bytesize) if (r % 2 == 0) ; r -= 1 ; end if ( lower <= r && r <= higher && r.prime? == true ) return r end end end end class Integer def cdiv(d) # ceil divide w=self.divmod(d) if w[1] > 0 ; w[0] +=1 ; end return w[0] end def powm(i,n) # square-and-multiply for exp modulo n if ( i == 0 ) return 1 end b=1 a=self if ( (i & 1) == 1 ) b=a end i = i>>1 while ( i > 0 ) a = (a * a) % n if ( (i & 1) == 1 ) b = (a * b) % n end i = i>>1 end return b end def gcd(b) # GCD a = self if ( a < b ); t=b; b=a; a=t; end while b > 0 r = a % b a = b b = r end return a end def invert(n) # Extension Euclid Algorithm a = n # See Knuth's The Art of Computer Programming b = self % n # Vol2 pp.342 -- 343 p = 0 ; q = 1; v = 0; u = 1; while q > 0 p = a / b q = a % b w = v - (p*u) ## DEBUG printf "%d = %d(%d) + %d (%d)\n", a,p,b,q,v ### a = b b = q v = u u = w end return (v+n) % n end # Miller-Rabin Test (Prime Test) # See, http://www.cs.albany.edu/~berg/csi445/Assignments/pa4.html def bitarray(n) b=Array::new i=0 v=n while v > 0 b[i] = (0x1 & v) v = v/2 i = i + 1 end return b end def miller_rabin(n,s) b=bitarray(n-1) i=b.size j =1 while j <= s a = 1 + (rand(n-2).to_i) if witness(a,n,i,b) == true return false end j+=1 end return true end def witness(a,n,i,b) d=1 x=0 while i > 0 x = d d = (d**2) % n if ( (d == 1) && (x != 1) && (x != (n-1)) ) return true end if ( b[i-1] == 1 ) d = (d * a ) % n end i -= 1 end if ( d != 1) return true end return false end def prime? a = self return miller_rabin(a,30) end def bytesize n = self i=0 while n > 0 n = n >> 8 i += 1 end return i end def bitsize n = self i=0 while n > 0 n = n >> 1 i += 1 end return i end end