wheel factorizationとか言うやつ

wheel.rb

# M_kを、k番目までの素数の積を合計したものとする。
# 例えばM2ならば2 * 3で6となる。これがホイールの
# 周期として定義される。

# M2の場合、最初の素数以降、確率的に2, 4, 2, 4...
# の間隔で、素数の候補が表われる。これらは、自分の理解
# するところによれば、2と3の倍数が表れない周期性を利用
# している。

M2_CYCLE_LENGTH = 6
M2_CYCLE = [4, 2]

# 実際に生成されたリストを見ると
#
#  1,  5,
#  7, 11,
# 13, 17,
# 19, 23,
# 25, 29..
#
# となるわけだが、25は素数ではない。
#
# wheel factorizationのおそらくの弱点は、
# 任意の倍数から除外される合成数を見逃すところにある。
# 実際、25は2と3では割りきれない数である

# この解決方法は非常に簡単である。
# それは、この周期を拡大すればよい。M3 = 2 * 3 * 5 = 30
# に拡張する。

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

# 当然、M3の場合にも合成数があり、49 = 7 * 7 が該当する。

def wheel_factorization limit
  wheel_numbers = []
  next_number = 11

  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

# ---------------------------------------------
# 以下 Helper

def not_prime? x, y
  !(x == 1 || y == 1 || x == y) && (x % y == 0)
end

def print_not_prime numbers
  numbers.each { |x|
    numbers.each {|y| puts "#{x} = #{y} * #{x / y}" if not_prime? x, y}
  }
end

LIMIT = 10000
numbers = wheel_factorization LIMIT

# -------------------------------
# 素数だけにする
# real    0m5.688s
# user    0m5.672s
# sys     0m0.012s

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

primes = prime_check! numbers

#-------------------------------
# M3を手で求めるのはさすがにしんどい
def print_differences numbers
  prev = 0
  diffs = []
  numbers.each do |x|
    puts "#{x} = #{prev} + #{x - prev} "
    diffs << (x - prev)
    prev = x
  end
  diffs
end

# テスト用エラトステネスのふるい
def eratos limit
  numbers = [*2..limit]
  primes  = []
  until numbers == []
    number = numbers[0]
    numbers = numbers.select { |x| !(x % number == 0)}.to_a
    primes << number
  end
  primes
end
p eratos(LIMIT) - primes