ruby-dev

Yuguiです。


keiju ISHITSUKA さんは書きました:
> 何カ所か変更個所がありましたが, いかがでしょう? 
> 
> もし, 問題ないようであれば, すいませんが, それを変更してcomit していた
> だけないです?

ありがとうございます。コミットしようかと思いましたがテストを書いていたら
2点問題に気づきました。結果としてちょっと変更が必要になりそうですので、
恐れ入りますが再度確認をお願いします。

先にお送りしたprime.rbでは2点問題がありました。
1. Prime.eachが返すenumerator(実際にはgenerator)が、uboundを記憶しない。
   Prime.each(100){|p| }は停止するのに、Prime.each(100).each{|p| }は停止
しないという奇妙な結果になります。

2. Prime#nextがない。
  互換性にこだわって見せた割にうっかりPrime#nextの実装を書き忘れていました。
  また、先のバージョンを元に単純にPrime#nextを実装するとPrime#nextによっ
て進んだ内部ポインタとPrime#eachが同期しません。Prime#eachは先立つ
Prime#nextに関わらず2から列挙を開始してしまいます。
  一方、Prime::eachやPrime::instance.eachのほうは先立つPrimeのインスタン
スメソッド呼び出しには影響されずに2から列挙することを期待されると思います。


この2点の問題は解決できましたが、generatorのライフサイクルが少々変わりま
した。新しいprime.rbを添付します。これで差し支えないようでしたらPrimeを
削除したバージョンのmathn.rbと併せてコミットしようと思います。いかがで
しょうか?


-- 
Yugui <yugui@yugui.jp>
http://yugui.jp
私は私をDumpする

Attachments (2)

prime.rb (9.48 KB, text/x-ruby)

require "singleton"
require "forwardable"

class Integer
  # Re-composes a prime factorization and returns the product.
  #
  # See Prime#int_from_prime_division for more details.
  def Integer.from_prime_division(pd)
    Prime.int_from_prime_division(pd)
  end
  
  # Returns the factorization of +self+.
  # 
  # See Prime#prime_division for more details.
  def prime_division(generator = Prime::Generator23.new)
    Prime.prime_division(self, generator)
  end

  # Returns true if +self+ is a prime number, false for a composite.
  def prime?
    Prime.prime?(self)
  end

  # Iterates over all prime numbers. 
  #
  # Calls +block+ once for each prime numer, passing a prime as
  # a parameter.
  # +ubound+:: Upper bound of prime numbers. The iterator stops 
  #            after yields all prime numbers +p <= ubound+.
  def Integer.each_prime(ubound, &block) # :yields: prime
    Prime.each(ubound, &block)
  end
end

# The set of all prime numbers.
#
# == Example
#  Prime.each(100) do |prime|
#    p prime  #=> 2, 3, 5, 7, 11, ...., 97
#  end
class Prime
  include Enumerable

  def initialize # :nodoc:
    @generator = nil
  end
  @the_instance = Prime.new

  # obsolete.
  #
  # Use +Prime::instance+ or class methods of +Prime+.
  def initialize
    @generator = EratosthenesGenerator.new
    warn "Prime::new is obsolete. use Prime::instance or class methods of Prime."
  end

  def succ
    @generator.succ
  end
  alias next succ

  class<<self
    extend Forwardable
    include Enumerable
    # Returns the default instance of Prime.
    def instance; @the_instance end

    def method_added(method) # :nodoc:
      (class<<self;self;end).def_delegator :instance, method
    end
  end

  # iterates the given block over all prime numbers.
  def each(ubound = nil, generator = nil, &block)
    generator ||= (@generator || EratosthenesGenerator.new)
    orig_ubound = generator.upper_bound
    generator.upper_bound = ubound
    generator.each(&block)
  ensure
    generator.upper_bound = orig_ubound
  end


  # returns true if +value+ is prime, false for a composite.
  #
  # +value+:: an arbitrary integer.
  # +generator+:: optional. A pseudo-prime generator.
  def prime?(value, generator = Prime::Generator23.new)
    for num in generator
      return true if value < num * num
      return false if value % num == 0
    end
  end

  # Re-composes a prime factorization and returns the product.
  #
  # +pd+:: Array of arrays of integers. The each internal 
  #        array consists of a prime number and a natural number. 
  #
  # For +[[p_1, e_1], [p_2, e_2], ...., [p_n, e_n]]+, it returns
  # +p_1**e_1 * p_2**e_2 * .... * p_n**e_n+.
  def int_from_prime_division(pd)
    pd.inject(1){|value, (prime, index)|
      value *= prime**index
    }
  end

  # Returns the factorization of +value+.
  #
  # For an arbitrary integer 
  # +n = p_1**e_1 * p_2**e_2 * .... * p_n**e_n+,
  # +prime_division(n)+ returns
  # +[[p_1, e_1], [p_2, e_2], ...., [p_n, e_n]]+.
  #
  # === Parameters
  # +value+:: An arbitrary integer.
  # +generator+:: Optional. A pseudo-prime generator.
  #               +generator.succ+ must return the next 
  #               pseudo-prime number in the ascendent
  #               order. It must generate all prime numbers,
  #               but may generate non prime numbers.
  #
  # ==== Exceptions
  # +ZeroDivisionError+:: when +value+ is zero.
  def prime_division(value, generator= Prime::Generator23.new)
    raise ZeroDivisionError if value == 0
    pv = []
    for prime in generator
      count = 0
      while (value1, mod = value.divmod(prime)
	     mod) == 0
	value = value1
	count += 1
      end
      if count != 0
	pv.push [prime, count]
      end
      break if prime * prime  >= value
    end
    if value > 1
      pv.push [value, 1]
    end
    return pv
  end

  # An abstract class for enumerating pseudo-prime numbers.
  #
  # Concrete subclasses should override succ, next, rewind.
  class PseudoPrimeGenerator
    include Enumerable

    def initialize(ubound = nil)
      @ubound = ubound
    end

    def upper_bound=(ubound)
      @ubound = ubound
    end
    def upper_bound
      @ubound
    end

    # returns the next pseudo-prime number, and move the internal
    # position forward. 
    #
    # +PseudoPrimeGenerator#succ+ raises +NotImplementedError+. 
    def succ
      raise NotImplementedError, "need to define `succ'"
    end

    # alias of +succ+.
    def next
      raise NotImplementedError, "need to define `next'"
    end

    # Rewinds the internal position for enumeration.
    #
    # See +Enumerator#rewind+.
    def rewind
      raise NotImplementedError, "need to define `rewind'"
    end

    # Iterates the given block for each prime numbers.
    # +ubound+:: 
    def each(&block)
      return self.dup unless block
      if @ubound
	loop do
	  p = succ
	  break if p > @ubound
	  block.call p
	end
      else
	loop do
	  block.call succ
	end
      end
    end

    # see +Enumerator#with_index+.
    alias with_index each_with_index

    # see +Enumerator#with_object+.
    def with_object(obj)
      return enum_for(:with_object) unless block_given?
      each do |prime|
	yield prime, obj
      end
    end
  end
  
  # An implementation of +PseudoPrimeGenerator+.
  #
  # Uses +EratosthenesSieve+.
  class EratosthenesGenerator < PseudoPrimeGenerator
    def initialize
      @last_prime = nil
    end
    
    def succ
      @last_prime = @last_prime ? EratosthenesSieve.instance.next_to(@last_prime) : 2
    end
    def rewind
      initialize
    end
    alias next succ
  end

  # An implementation of +PseudoPrimeGenerator+ which uses 
  # a prime table generated by trial division.
  class TrialDivisionGenerator<PseudoPrimeGenerator
    def initialize
      @index = -1
    end
    
    def succ
      TrialDivision.instance[@index += 1]
    end
    def rewind
      initialize
    end
    alias next succ
  end

  # Generates all integer which are greater than 2 and
  # are not divided by 2 nor 3.
  #
  # This is a pseudo-prime generator, suitable on 
  # checking primality of a integer by brute force 
  # method.
  class Generator23<PseudoPrimeGenerator
    def initialize
      @prime = 1
      @step = nil
    end
    
    def succ
      loop do
	if (@step)
	  @prime += @step
	  @step = 6 - @step
	else
	  case @prime
	  when 1; @prime = 2
	  when 2; @prime = 3
	  when 3; @prime = 5; @step = 2
	  end
	end
	return @prime
      end
    end
    alias next succ
    def rewind
      initialize
    end
  end




  # An implementation of prime table by trial division method.
  class TrialDivision
    include Singleton

    def initialize # :nodoc:
      # These are included as class variables to cache them for later uses.  If memory
      #   usage is a problem, they can be put in Prime#initialize as instance variables.

      # There must be no primes between @primes[-1] and @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]
      # @next_to_check % 6 must be 1.  
      @next_to_check = 103            # @primes[-1] - @primes[-1] % 6 + 7
      @ulticheck_index = 3            # @primes.index(@primes.reverse.find {|n|
      #   n < Math.sqrt(@@next_to_check) })
      @ulticheck_next_squared = 121   # @primes[@ulticheck_index + 1] ** 2
    end

    # Returns the cached prime numbers.
    def cache
      return @primes
    end
    alias primes cache
    alias primes_so_far cache

    # Returns the +index+th prime number. 
    #
    # +index+ is a 0-based index.
    def [](index)
      while index >= @primes.length
	# Only check for prime factors up to the square root of the potential primes,
	#   but without the performance hit of an actual square root calculation.
	if @next_to_check + 4 > @ulticheck_next_squared
	  @ulticheck_index += 1
	  @ulticheck_next_squared = @primes.at(@ulticheck_index + 1) ** 2
	end
	# Only check numbers congruent to one and five, modulo six. All others

	#   are divisible by two or three.  This also allows us to skip checking against
	#   two and three.
	@primes.push @next_to_check if @primes[2..@ulticheck_index].find {|prime| @next_to_check % prime == 0 }.nil?
	@next_to_check += 4
	@primes.push @next_to_check if @primes[2..@ulticheck_index].find {|prime| @next_to_check % prime == 0 }.nil?
	@next_to_check += 2 
      end
      return @primes[index]
    end
  end

  # An implementation of eratosthenes's sieve
  class EratosthenesSieve
    include Singleton

    def initialize # :nodoc:
      # bitmap for odd prime numbers less than 256.
      # For an arbitrary odd number n, @table[i][j] is 1 when n is prime where i,j = n.divmod(32) .
      @table = [0xcb6e, 0x64b4, 0x129a, 0x816d, 0x4c32, 0x864a, 0x820d, 0x2196]
    end

    # returns the least odd prime number which is greater than +n+.
    def next_to(n)
      n = (n-1).div(2)*2+3 # the next odd number of given n
      i,j = n.divmod(32)
      loop do
	extend_table until @table.length > i
	if !@table[i].zero?
	  (j...32).step(2) do |j|
	    return 32*i+j if !@table[i][j.div(2)].zero?
	  end
	end
	i += 1; j = 1
      end
    end

    private
    def extend_table
      orig_len = @table.length
      new_len = [orig_len**2, orig_len+256].min
      lbound = orig_len*32
      ubound = new_len*32
      @table.fill(0xFFFF, orig_len...new_len)
      (3..Integer(Math.sqrt(ubound))).step(2) do |p|
	i, j = p.divmod(32)
	next if @table[i][j.div(2)].zero?

	start = (lbound.div(2*p)*2+1)*p    # odd multiple of p which is greater than or equal to lbound
	(start...ubound).step(2*p) do |n|
	  i, j = n.divmod(32)
	  @table[i] &= 0xFFFF ^ (1<<(j.div(2)))
	end
      end
    end
  end
end

test_prime.rb (2.72 KB, text/x-ruby)

require 'test/unit'
require 'prime'
require 'stringio'

class TestPrime < Test::Unit::TestCase
  # The first 100 prime numbers
  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, 241, 251, 257, 263, 269, 271, 277, 
    281, 283, 293, 307, 311, 313, 317, 331, 337, 
    347, 349, 353, 359, 367, 373, 379, 383, 389, 
    397, 401, 409, 419, 421, 431, 433, 439, 443, 
    449, 457, 461, 463, 467, 479, 487, 491, 499, 
    503, 509, 521, 523, 541,
  ]
  def test_each
    primes = []
    Prime.each do |p|
      break if p > 541
      primes << p
    end
    assert_equal PRIMES, primes
  end

  def test_each_by_prime_number_theorem
    3.upto(15) do |i|
      max = 2**i
      primes = []
      Prime.each do |p|
        break if p >= max
        primes << p
      end

      # Prime number theorem
      assert primes.length >= max/Math.log(max)  
      assert primes.length <= (2..max).step(0.1).inject(0){|sum,x| sum + 0.1 / Math.log(x)} # Li(x)

#      assert Prime.cache.length >= primes.length
    end
  end

  def test_each_without_block
    enum = Prime.each
    assert enum.respond_to?(:each)
    assert enum.kind_of?(Enumerable)
    assert enum.respond_to?(:with_index)
    assert enum.respond_to?(:next)
    assert enum.respond_to?(:succ)
#    assert enum.respond_to?(:rewind)
  end

  def test_new
    buf = StringIO.new('', 'w')
    orig, $stderr = $stderr, buf

    enum = Prime.new
    assert !buf.string.empty?
    $stderr = orig

    assert enum.respond_to?(:each)
    assert enum.kind_of?(Enumerable)
    assert enum.respond_to?(:succ)

    assert Prime === enum
  ensure
    $stderr = orig
  end

  def test_enumerator_succ
    enum = Prime.each
    assert_equal PRIMES[0, 50], 50.times.map{ enum.succ }
    assert_equal PRIMES[50, 50], 50.times.map{ enum.succ }
    enum.rewind
    assert_equal PRIMES[0, 100], 100.times.map{ enum.succ }
  end

  def test_enumerator_with_index
    enum = Prime.each
    last = -1
    enum.with_index do |p,i|
      break if i >= 100
      assert_equal last+1, i
      assert_equal PRIMES[i], p
      last = i
    end
  end

  class TestInteger < Test::Unit::TestCase
    def test_prime_division
      pd = PRIMES.inject(&:*).prime_division
      assert_equal PRIMES.map{|p| [p, 1]}, pd
    end

    def test_from_prime_division
      assert_equal PRIMES.inject(&:*), Integer.from_prime_division(PRIMES.map{|p| [p,1]})
    end

    def test_prime?
      assert 2.prime?
      assert 3.prime?
      assert (2**31-1).prime?
      assert !(2**32-1).prime?
      assert !((2**13-1) * (2**17-1)).prime?
    end
  end
end

Thread

Prev Next

In This Thread

Prev Next