読者です 読者をやめる 読者になる 読者になる

Line 1: Error: Invalid Blog('by Esehara' )

または私は如何にして心配するのを止めてバグを愛するようになったか

>> Zanmemo

あと何かあれは 「esehara あっと じーめーる」 か @esehara まで

Rubyで、「エラトステネスのふるい」を改良したと言われる「アトキンのふるい」をいまさら実装する

今日の料理

f:id:nisemono_san:20160410191310j:plain

ニラを入れすぎた

概要

素数のリストを作成する際に、エラトステネスのふるいが有名だが、それを改良した「アトキンのふるい」というものがある。今回は、他のコードを参考にしながら、自前で実装してみて、思ったことを書き残しておきたいと思う。

はじめに

前に素因数分解のコードを練習がてら書いた。そのときに、素数を列挙するために「エラトステネスのふるい」を利用した方法を使っていた。

「エラトステネスのふるい」は、非常に直感的かつ軽量なアプローチであると思う。なぜなら、素数の定義が「それ自身と1だけが約数になるような数」であるとするならば、「一度でもなんらかの数で割りきれれば、それは素数ではない」ということになる。ならば、「愚直に出てくる数で割っていけば良い」ということになる。

f:id:nisemono_san:20160414191033g:plain

日本のWikipidiaより参照

とはいえ、さすがに「この単純なアプローチ」が素数のリストを求めるのにベターであるとは言えないだろう。もう少し、このようなリストを探すのに、いい方法がないか探してみた結果、「アトキンのふるい」という方法があるということを、このブログで知った

ちなみに、こちらの方では改良版もあるので、この記事を読むよりいいと思う(自分の記事を先に否定する奴)。

概説

詳しくはこの論文が発端になっているようだ。この論文によれば、素数の出てくる周期性みたいなものを、三つの区間に分類していると考えるといいようだ。

一つ目の式

まず最初の区間として、{p \in 1 + 4Z \iff 4x^{2} + y^{2} = n}の正の解(x, y)が奇数個存在していることが条件になる。このとき、正の解を何如にして求めるか、ということになる。

まず、{ p \in 1 + 4Z }ということに関しては、他のコードを見る限り、「12で割って、1の余りか、5の余りになるもの」という実装をされていることが多い。(ちなみに、「4で割って1が余るもの」として厳密に考えた場合、9も上記に該当するが、これは素数ではないとわかるので排除されているようだ)

そして、問題は{4x^{2} + y^{2} = n}の解である(x, y)を求めることだが、これに関しては、ループを廻して総当たりにするようだ。

class SieveHelper
  # .....
   def first_formula
    (4 * @x ** 2) + (@y ** 2)
  end

  def first_check?
    n = first_formula
    n <= @limit && (n % 12 == 1 || n % 12 == 5)
  end

  def puts_first_formula
    puts "#{first_formula}, x = #{@x}, y = #{@y}"
  end
end

このようなヘルパーを用意し、

limit = 100
sqrt_limit = Math.sqrt(limit).ceil
sieve = Array.new(limit)

1.upto(sqrt_limit) do |x|
  1.upto(sqrt_limit) do |y|
    h = SieveHelper.new(limit, x, y)
    if h.first_check?
      n = h.first_formula
      h.puts_first_formula
      sieve[n] = !sieve[n]
    end
  end
end

といったようなループを回すことにする。奇数個あったかどうかをわざわざ数えなくても、奇数というのは、2間隔毎にあらわれるわけだから、最初の0回をfalseとし、あとはこれを交互に切りかえることで、数えられることができる。

以下、100までの素数候補のリストである(2回出てきたものについては省略した)。

5, x = 1, y = 1
13, x = 1, y = 3
17, x = 2, y = 1
25, x = 2, y = 3
29, x = 1, y = 5
37, x = 3, y = 1
41, x = 2, y = 5
53, x = 1, y = 7
61, x = 3, y = 5
73, x = 4, y = 3
89, x = 4, y = 5
97, x = 2, y = 9

ちなみに、「25」という、あきらかに5の素数のものが入っているが、これについては、のちに、これで候補を絞りこんだ数のリストで二重にチェックすることになる。

また、素数リストとして、これを眺めた場合、気がつくのは、「5, 13, 17, 29, 37...」の間に、「7, 11, 19, 23, 31」が無いということになる。なので、2番目の式で補足することにする。

2つ目の式

1つ目の式で、だいたいの実装の方針はわかったと思う。次の式は{p \in 7 + 12Z \iff 3x^{2} + y^{2} = n}の正の解(x, y)が奇数個存在していることが条件になる。なので、今回も素朴に実装する。

  #...
  
  def second_formula
    (3 * @x ** 2) + (@y ** 2)
  end

  def second_check?
    n = second_formula
    n <= @limit && (n % 12 == 7)
  end

さて、この場合だと次のような解が得られる。

7, x = 1, y = 2
19, x = 1, y = 4
31, x = 3, y = 2
43, x = 3, y = 4
67, x = 1, y = 8
79, x = 5, y = 2

さて、先ほどの「5, 13, 17, 29, 37...」のリストでは、「7, 11, 19, 23, 31」が足りなかった、これによって「5, 7, 13, 19, 17, 29, 31, 37...」が揃ったということになる。

さて、見る限り、まだ「11, 23」が出てきていない。こいつも補足しなくてはならない。

3つ目の式

さっそく実装を行なう。3つ目は{p \in 11 + 12Z \iff 3x^{2} - y^{2} = n, x > y > 0}の正の解(x, y)が奇数個存在していることが条件になる。

  # ...

  def third_formula
    (3 * @x ** 2) - (@y ** 2)
  end

  def third_check?
    n = third_formula
    n <= @limit && @x > @y && (n % 12 == 11)
  end

さて、これらの結果は以下の通りである。

11, x = 2, y = 1
23, x = 3, y = 2
47, x = 4, y = 1
71, x = 5, y = 2
59, x = 5, y = 4
83, x = 6, y = 5

この段階で、「11, 23」が出てくる。これによって、少なくとも、30以下の整数に関しては、素数が補足できることが確認された。

しかし、とはいえもう少し問題がある。というのは、最初の式において「25」という「素数ではない数」が紛れこんでいた。ということは、他の式でも「25」のような、余計な数が入りこんでいないとはいい切れない。従って、これを排除するような関数を定義する必要がある。

最後のふるいわけ

def delete_if_multiple(sieve, limit)
  for n in 5..limit do
    if (sieve[n])
      x = n ** 2
      x.step(limit, x) {|i| sieve[i] = false }
    end
  end
  return sieve
end

さきほと、残っていたのが「25」であることを確認した。上の関数は、それぞれの数の2乗を消していくものである。これによって、当然「25」は排除される。

まとめ

これが自分が理解した範囲でのアトキンのふるいの実装 である。証明に関しては、なかなか複雑のようで、自分には手に負えないものであった。識者がいれば、解説を書いて欲しいと思う。

また、このコードについてはこのgistに掲載したので、もし参考にする奇特な人があれば、参考にして欲しい。

最後にリファレンスとして、参考にしたコードを掲載する。

こういう地味ながら知見が詰まれたアルゴリズムを実装するのは、とても勉強になるので、今後もやっていきたい。