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

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

>> Zanmemo

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

素数候補を効率よく絞りこむ「Wheel factorization」について調べたこと

今日の料理

f:id:nisemono_san:20160415004114j:plain

サーモン漬け丼

概要

素数リストの求め方に関しては、前回アトキンのふるいを紹介した。この他にもWheel factorizationという方法がある。調べてみたところ、直接的な日本語の紹介が乏しいので、実装したついでに、調べた範囲のことについてメモをしておく。

はじめに

素数リストを求める方法にも様々あって、前回はアトキンのふるいを紹介したのだが、論文を読めばわかるように、確かに高速なのかもしれないが、あまり直感的とは言いがたい側面がある。

そんな中、アトキンのふるいに関して調べていた最中にJulia の素数周りの実装 - Qiitaという記事を見つけた。これによれば、「Wheel factorizationを使った Eratosthenesのふるい」にアルゴリズムを変更したことが書いてある。しかし、このWheel factorizationの日本語情報を探してみたのだが、自分の探しかたが悪かったのだろうか、殆ど見つからなかった。なので、必死にこういった解説を読みながら、理解につとめていたりしていた。

Wheel factorizationの発想について

f:id:nisemono_san:20160415144627p:plain

英語版Wikipediaより。

Wheel factorizationは、エラトステネスのふるいと合わせて使われる。何故ならば、Wheel factorization自体は、単純に「このあたりに素数がありそうだ」という絞りこみを行なうためのものだからである。しかし、この絞りこみは、総当たりで試し割りするよりも計算量を抑えることができる。

その絞りこみの方法だけれども、素数があらわれそうな周期性に注目しているといってもいい。ここで、「ある順番までの素数の積」を{M_i}として定義する。話を簡単にするために、{M_2}とする。このとき、2番目までの素数の積は2 * 3 = 6となる。

この2 * 3 = 6が一つのサイクルとなる。このとき、6つ毎に列を作ってみる。

 1, 2, 3, 4, 5, 6,
 7, 8, 9,10,11,12,
13,14,15,16,17,18..

まず、サイクル最後の数は明らかに先程素数同士をかけあわせた数なので、素数ではない、なのでこれらは排除できる。次に、このサイクルを見ると、2, 3の縦の列に関して倍数であることがわかる。従ってこれも除去ができる。当然、4も同様に2の倍数なので、これを除去しよう。

 1, 2, 3,  , 5,  ,
 7,  ,  ,  ,11,  ,
13,  ,  ,  ,17,  ..

さて、ついでに最初に使った素数も除去してみると、不思議なことが起きる。

 1,  ,  ,  , 5,  ,
 7,  ,  ,  ,11,  ,
13,  ,  ,  ,17,  ..

このように、素数が表れるサイクルが見つかる!

とはいえ、これは最初だから上手くいくのであって、表を続けてみればわかるのだけれども、素数ではない数があらわれる。5列目の最初の数は25であり、これは5の2乗なので、素数ではない。

このサイクルの弱点は、最初に使った素数以降に見付かった素数で割れるような合成数を排除できないことである。今回ならば25である。とはいえ、調べる範囲が3分の1になっただけでも十分である。

とすると、このサイクルを大きく取れば、もっと効率よくなるのではないか、と考えるのは常である。そこで、実際にWheel factorizationを実装してみることにする。

実装

まず最初に{M_3}のデータを与える。

M3_CYCLE_LENGTH = 30
M3_PRIMES = [2, 3, 5]
M3_CYCLE = [6, 4, 2, 4, 2, 4, 6, 2]

この周期に従って、素数候補リストを生成する。

def wheel_factorization limit
  wheel_numbers = []
  next_number = 1

  1.upto(limit) do |n|
    wheel_numbers << next_number
    next_number += M3_CYCLE[n % M3_CYCLE.size - 1]
    break if limit < next_number
  end
  return M3_PRIMES + wheel_numbers - [1]
end

あとは、エラトステネスのふるいを使って素朴に素数っでないものを排除すればよい。

def prime_check! numbers
  primes = []
  until numbers == []
    number = numbers[0]
    numbers = numbers.select {|x| !(x % number == 0)}.to_a
    primes << number
  end
  primes
end

議論

一見、このサイクルを大きくすればするほど、効率的に見えるかもしれない。しかし、2 * 3 * 5 * 7 = 210になると不幸なことが起きる。

ふるいを使うということは、「素数でないものが余計に入っている場合」はいいのだけれども、「素数だったものがふるい落とされてしまう」とどうしようもない。

実際に、2-3-5-7にしてしまうと、331を見落としてしまうという問題が指摘されている。ちなみに先程のJulia langにおいても30サイクルを利用しているようで、大きくすればいいというわけではないようだ。

他のリファレンス

まとめ

人類が滅びて猫が文明を持ったときにも素数に興味を持つのか、ということに関心が出てきてしまうので困っていつつある。