Prime::A2EGenerator < Prime::PseudoPrimeGenerator

prime.a2e.rb

# original: Prime::EratosthenesGenerator and Prime::EratosthenesSieve
require 'prime'
require 'singleton'

class Prime

  # An implementation of +Prime::PseudoPrimeGenerator+.
  #
  # Uses +A2ESieve+.
  class A2EGenerator < PseudoPrimeGenerator
    def initialize
      @last_prime = nil
      super
    end

    def succ
      case @last_prime
        when nil
          @last_prime = 2
        when 2
          @last_prime = 3
        else
          @last_prime = A2ESieve.instance.next_to(@last_prime)
      end
    end
    alias next succ

    def rewind
      initialize
    end
  end

  # Internal use. An implementation of eratosthenes's sieve
  class A2ESieve
    include Singleton

    def initialize # :nodoc:
      # bitmap for odd prime numbers less than 384.
      # For an arbitrary odd number n except for the multiples of 3, 
      # @table[i][j] is 1 when n is prime where i,j = n.div(3).divmod(16).
      @table = [0xf6fe, 0x2dda, 0x6c3f, 0x9ad6, 0xbc47, 0x76a9, 0x43c2, 0xd4b9]
    end

    # returns the least odd prime number which is greater than +n+.
    def next_to(n)
      i, j = ((n-1).div(3)+1).divmod(16)
      while true
        extend_table until @table.length > i
        if !@table[i].zero?
          j.upto(15) do |k|
            return i*48 + k*3 + k%2 + 1 if !@table[i][k].zero?
          end
        end
        i += 1; j = 0
      end
    end

    private
    def extend_table
      orig_len = @table.length
      new_len = [orig_len**2, orig_len+256].min
      lbound = orig_len * 48
      ubound = new_len * 48
      @table.fill(0xffff, orig_len...new_len)
      p = 1
      p_ubound = Integer(Math.sqrt(ubound))
      d = 4
      while (p += d) <= p_ubound
        d = 6 - d
        i, j = p.div(3).divmod(16)
        next if @table[i][j].zero?
        # odd multiple of p and not of 3 which is nearly equal to lbound
        n = ((lbound + 5*p).div(6*p) * 6 - 1) * p
        s = 2
        while n < ubound
          i, j = n.div(3).divmod(16)
          @table[i] &= 0xffff ^ (1<<j)
          n += s * p
          s = 6 - s
        end
      end
    end
  end
end

if $0 == __FILE__
  max = ARGV.size == 1 ? ARGV[0].to_i : 100
  pr = 0
  Prime.each(max, Prime::A2EGenerator.new){|num|pr=num}
  puts "The max prime number (<=#{max}) is #{pr}."
  pr = Prime.each(nil, Prime::A2EGenerator.new).take(max)[-1]
  puts "The #{max}th prime number is #{pr}."
end