Project Euler 200 を真面目に解く

以前、「2, 5 の倍数以外は、耐素数性はない」と仮定して解いたが、200番目までは問題ないようだが、ずっとそうではないようなので、完全に解けるアルゴリズムを考えた。

  • まず a**2, b**3 (a,bは素数)のそれぞれの1次配列を作成する
    • lista = [4, 9, 25, 49, 121, ...]
    • listb = [8, 27, 125, 343, 1331, ...]
  • ループ回数は a, b の全パターン。ただし、a**2 * b**3 の小さい順にループする。すなわち、200番目の条件を満たすsqubeが見つかるまで。(ちなみに、これは 43898 ループとなる)
    • ここで、「a**2 * b**3 の小さい順にループする」というのがミソ
  • ループ中に、listb の各項に対して、lista のどの項まで検査したかを保存するための1次配列(jdglist)も用意する
    • 例えば、a**2 * b**3 = 1125 まで検査したときの状況は以下となる
      • listb[0] は lista[3] まで検査済み
      • listb[1] は lista[2] まで検査済み
      • listb[2] は lista[1] まで検査済み
      • listb[3] は lista[0] まで検査済み
      • listb[4] は lista[0] まで検査済み
    • jdglist の中は以下のようになる
      • jdglist = [3, 2, 1, 0, 0, nil, ...]
    • この状態で、次に検査する sqube は何になるか?は、そのループ中に判断できて、以下のうち一番小さい数字となる (※1)
      • listb[0] * lista[4]
      • listb[1] * lista[3]
      • listb[2] * lista[2]
      • listb[3] * lista[1]
      • listb[4] * lista[1]
      • listb[5] * lista[0]


このアルゴリズムで計算すると、200番目を求める時間は 56s となる。
ちなみに、面白いことに、jdglist を「lista の各項に対して、listb のどの項まで検査したかを保存するための1次配列」と定義すると、835s かかる。jdglist のサイズが大きくなり、(※1) の計算時間が何倍も大きくなるからだと思われる。


ということで、ループの中はこんな感じになりました。

  loop do
    ### a**2 * b**3 の小さい数字から検査する

    if listb.size == jdglist.size
      puts "over listb. #{jdglist.size}"
      return
    end

    templist = []

    (0..jdglist.size-1).each do |itcur|
      b = lista[jdglist[itcur]+1]
      if b.nil?
        puts "over lista. #{jdglist[itcur]+1}"
        return
      end
      templist << [listb[itcur] * b, itcur]
    end

    templist << [listb[jdglist.size] * lista[0], jdglist.size]

    min = templist.min_by {|t| t[0]}

    if jdglist[min[1]].nil?
      jdglist[min[1]] = 0
    else
      jdglist[min[1]] += 1
    end

    ### 条件に合うものを selected に入れる

    sqube = min[0]

    if include200(sqube) && prime_proof.is?(sqube)
      count += 1
      puts "#{count}: (a, b)=(#{min[1]}, #{jdglist[min[1]]}) sqube = #{sqube}"
      selected << sqube
    end

    break if count == 200
  end