Project Euler 41、高速な素数判定
Pandigitalな数を作り、それを1つずつ素数判定する方法で求めました。(まあ、素直で愚直な方法です)
今まで素数判定について、エラトステネスの篩で求めた素数に含まれているか?という手抜き方法で実装してきましたが、流石に速度をごまかせなくなったので高速な素数判定アルゴリズムを調べることにしました。
ちなみにエラトステネスの篩で現実的に素数判定できるのは、Rubyでは100万までですね。1000万レベルは、たぶんかなりキツイ。1分近くかかってしまいます。今回は9桁なのでエラトステネスの篩では無理です。
ということで、調べた結果利用することにしたのが「ミラー・ラビン素数判定法」。
とはいうものの、結局アルゴリズムは理解できず、ソースコードはネットにあるものをそのまま流用した感じになってしまいました。。。
やっぱり数学は難しいでつ。 つД`)
ということで参考にさせていただいたサイトはこちら。
- http://d.hatena.ne.jp/pashango_p/20090704/1246692091
- http://www.weblio.jp/wkpja/content/%E3%83%9F%E3%83%A9%E3%83%BC-%E3%83%A9%E3%83%93%E3%83%B3%E7%B4%A0%E6%95%B0%E5%88%A4%E5%AE%9A%E6%B3%95_%E3%82%B3%E3%83%BC%E3%83%89%E4%BE%8B
Pythonには pow(a,b,c) という、a^b mod c を求めるメソッドが標準に用意されていて、これがすごく早い。おそらく「モンゴメリ乗算」がC言語で実装されているのではないかとのこと。
Ruby にはそのようなメソッドは標準では用意されておらず、自作せざるを得ず、結局「ミラー・ラビン素数判定法」はPythonの方が圧倒的に早いという結果になってしまいました。(目測3倍)
ということで今回用意した素数判定メソッドはこちら。
module ModMath def self.pow(base, power, mod) result = 1 while power > 0 result = (result * base) % mod if power & 1 == 1 base = (base * base) % mod power >>= 1 end result end def self.is_prime(q, k=50) q = q.abs return true if q == 2 return false if q < 2 || q&1 == 0 d = (q-1)>>1 while d&1 == 0 do d >>= 1 end k.times do a = rand(q-2) + 1 t = d y = ModMath.pow(a, t, q) while t != q-1 && y != 1 && y != q-1 do y = ModMath.pow(y, 2, q) t <<= 1 end return false if y != q-1 && t&1 == 0 end return true end end
そしてProblem41のコードはこれです。(ちょー適当です)
def problem41 def getmaxprime(aryall) maxprime = -1 firstlist = aryall.dup firstlist.delete_if {|t| t!=3 && t!=7 && t!=9} firstlist.each do |first| aryalldup = aryall.dup aryalldup.delete_if {|t| t == first} permutation(aryalldup).each do |pand| num = digits_to_num([first] + pand) if ModMath.is_prime(num) p num maxprime = (maxprime < num) ? num : maxprime end end end maxprime end maxprime = -1 maxprime = getmaxprime([1,2,3,4,5,6,7,8,9]) if maxprime < 0 maxprime = getmaxprime([1,2,3,4,5,6,7,8]) if maxprime < 0 maxprime = getmaxprime([1,2,3,4,5,6,7]) if maxprime < 0 maxprime = getmaxprime([1,2,3,4,5,6]) if maxprime < 0 maxprime = getmaxprime([1,2,3,4,5]) if maxprime < 0 maxprime = getmaxprime([1,2,3,4]) if maxprime < 0 maxprime end
あ、そういえば順列を作るメソッドを自作してたんだった。あと、各桁を並べた数列を数字に戻すメソッドも。
def permutation(ary) return [ary] if ary.size == 1 res = [] (0..ary.size-1).each do |t| arydup = ary.dup arydup.delete_at(t) permutation(arydup).each do |subary| res << [ary[t]] + subary end end res end def digits_to_num(dig) num = 0 t = 0 while dig.size > t num += dig[t] * 10**t t += 1 end return num end
9桁のPandigital数の素数判定に 16.8125s もかかっており、最終的にプログラムが終わるのも 19.25s と結構時間がかかる。
まあ、Project Euler の "one-minute rule" は満たしているからいいけど、たぶんもっと早いアルゴリズムはあるんだろうな。