Project Euler 41、高速な素数判定

プログラミング 数学

Pandigitalな数を作り、それを1つずつ素数判定する方法で求めました。（まあ、素直で愚直な方法です）

今まで素数判定について、エラトステネスの篩で求めた素数に含まれているか？という手抜き方法で実装してきましたが、流石に速度をごまかせなくなったので高速な素数判定アルゴリズムを調べることにしました。

ちなみにエラトステネスの篩で現実的に素数判定できるのは、Rubyでは100万までですね。1000万レベルは、たぶんかなりキツイ。1分近くかかってしまいます。今回は9桁なのでエラトステネスの篩では無理です。


ということで、調べた結果利用することにしたのが「ミラー・ラビン素数判定法」。
とはいうものの、結局アルゴリズムは理解できず、ソースコードはネットにあるものをそのまま流用した感じになってしまいました。。。
やっぱり数学は難しいでつ。 つД`)


ということで参考にさせていただいたサイトはこちら。


Pythonには pow(a,b,c) という、a^b mod c を求めるメソッドが標準に用意されていて、これがすごく早い。おそらく「モンゴメリ乗算」がC言語で実装されているのではないかとのこと。

Ruby にはそのようなメソッドは標準では用意されておらず、自作せざるを得ず、結局「ミラー・ラビン素数判定法」はPythonの方が圧倒的に早いという結果になってしまいました。（目測３倍）


ということで今回用意した素数判定メソッドはこちら。

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" は満たしているからいいけど、たぶんもっと早いアルゴリズムはあるんだろうな。