151212(3)

Ruby


Hilbert prime(1)

オンライン整数列大辞典の
A057948（http://oeis.org/A057948/list）
と比較し、答え合わせしてみる。

=begin
p(=4k + 1)をHilbert primeとし、2p, 3p, …をfalseにしていく。
ちなみに、
2p, 3p, 4p, 6p, …は
2p = 4(2k) + 2,
3p = 4(3k) + 3,
4p = 4(4k + 1),
6p = 4(6k + 1) + 2,
…
なので、listを用意した時点でfalseになっている。
=end

def A057948(n)
  # m以下のHilbert primeを探す
  m = 10 ** 6
  list = Array.new(m + 1, false)
  1.step(m, 4){|i| list[i] = true}
  i = 5
  while i * i <= m
    if list[i]
      j = i + i
      while j <= m
        list[j] = false
        j += i
      end
    end
    i += 1
  end
  (5..m).select{|i| list[i]}[0..n]
end
ary = A057948(55)

# OEIS A057948のデータ
ary0 =
[5,9,13,17,21,29,33,37,41,49,53,57,61,69,73,77,89,
 93,97,101,109,113,121,129,133,137,141,149,157,161,
 173,177,181,193,197,201,209,213,217,229,233,237,
 241,249,253,257,269,277,281,293,301,309,313,317,
 321,329]
# 一致の確認
p ary == ary0