From: C Erler <erlercw@...>
Date: 2006-09-28T12:21:47+09:00
Subject: Prime sieve optimization.

The prime sieve of Atkin (http://en.wikipedia.org/wiki/Sieve_of_Atkin)
is, if done right, faster than the sieve of Eratosthenes.  I've
included a reimplementation of the mathn Prime class using this sieve,
but it's slower than the one in ruby 1.9.

If anyone wants to try their hand at improving this (or starting anew
if this is a bad start), feel free.  Perhaps we can beat the one in
ruby 1.9 so we can replace it before 1.9 (1.10 ?) becomes the stable
version.

class Prime_Atkin
  include Enumerable

  # @@next_to_check is a multiple of 12.
  @@next_to_check = 240
  # @@primes should contain all primes in 1...@@next_to_check
  @@primes = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,
      73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,
      163,167,173,179,181,191,193,197,199,211,223,227,229,233,239]
  # @@primes_squared = @@primes.map { |prime| prime**2 }[2..-1]
  @@primes_squared = [25,49,121,169,289,361,529,841,961,1369,1681,
      1849,2209,2809,3481,3721,4489,5041,5329,6241,6889,7921,9409,
      10201,10609,11449,11881,12769,16129,17161,18769,19321,22201,
      22801,24649,26569,27889,29929,32041,32761,36481,37249,38809,
      39601,44521,49729,51529,52441,54289,57121]
  # SieveWidth is a multiple of 12.
  # Ensure that @@next_to_check + SieveWidth <= @@primes_squared.last
  SieveWidth = 56880
  @@new_primes = Array.new(SieveWidth, false)
  def initialize
    @index = -1
  end

  def succ
    @index += 1
    while @index >= @@primes.length
      range_end = @@next_to_check + SieveWidth
      high_x = Math.sqrt(range_end).floor

      1.upto(high_x) do |x|
        x_sq = x**2
        low_y = @@next_to_check - 4*x_sq
        if low_y < 1
          low_y = 1
        else
          low_y = Math.sqrt(low_y).ceil
        end
        high_y = [Math.sqrt(3*x_sq - @@next_to_check),
Math.sqrt(range_end - 3*x_sq)].max.floor

        high_y.downto(low_y) do |y|
          y_sq = y**2
          n = 4*x_sq + y_sq - @@next_to_check
          @@new_primes[n] = (not @@new_primes.at(n)) if (n >= 0 and n <
SieveWidth and (n % 12 == 1 or n % 12 == 5))
          n -= x_sq
          @@new_primes[n] = (not @@new_primes.at(n)) if (n >= 0 and n <
SieveWidth and n % 12 == 7)
          n -= 2*y_sq
          @@new_primes[n] = (not @@new_primes.at(n)) if (n >= 0 and n <
SieveWidth and x > y and n % 12 == 11)
        end
      end

      @@primes_squared.each do |prime_squared|
        multiple_index = prime_squared *
(@@next_to_check/prime_squared).ceil - @@next_to_check
        while multiple_index < SieveWidth
          @@new_primes[multiple_index] = false
          multiple_index += prime_squared
        end
      end

      @@new_primes.each_with_index do |prime_test, index|
        if prime_test
          prime = @@next_to_check + index
          @@primes << prime
          @@primes_squared << prime**2
        end
      end

      @@next_to_check = range_end
      @@new_primes.fill false
    end
    @@primes.at @index
  end
  alias next succ
  
  def each
    loop do
      yield succ
    end
  end
end