primes.rb

class Primes < Array
    def initialize(n)
        @max = n
        super n+1, false
        self[2] = true
        i = 3
        while i <= n
            self[i] = true
            i += 2
        end
        (3..n**0.5).each do |p|
            next unless self[p]
            unless self[p]
                p += 2
                next
            end
            q = p
            m = p*q
            while m <= n
                self[m] = false
                q += 2
                m = p*q
            end
        end
    end
    def seq
		Enumerator.new do |y|
			y << 2
			n = 3
			while n <= @max
				y << n if self[n]
				n += 2
			end
		end
    end
end

Prime_seq = Enumerator.new do |y|
    y << 2
    n = 3
    loop do
        ok = true
        p = 3
        while p*p <= n
            if n%p == 0
                ok = false
                break
            end
            p += 2
        end
        y << n if ok
        n += 2
    end
end

class Prime_factors < Hash
    def initialize(n)
        while n.even?
            self[2] ? self[2] += 1 : self[2] = 1
            n /= 2
        end
        p = 3
        while p*p <= n
            while n%p == 0
                self[p] ? self[p] += 1 : self[p] = 1
                n /= p
            end
            p += 2
        end
        self[n] = 1 if n > 1
    end
    def count
        self.values.inject(1) {|x,y| x*(y+1)}
    end
end

class Integer
    def is_prime      # "http://ja.wikipedia.org/wiki/ミラー-ラビン素数判定法" より、改変。
        n = self.abs()
        return true if n == 2
        return false if n == 1 || n.even?
        if n < 4_759_123_141 then
            arr = [2, 7, 61].take_while{|x| x < n}
        elsif n < 341_550_071_728_321 then
            arr = [2, 3, 5, 7, 11, 13, 17]
        else
            arr = []
            20.times do arr.push(rand(n-2)+1) end
        end
        d = n-1
        d /= 2 while d.even?
        arr.each do |a|
            t = d
            y = ModMath.pow(a,t,n)
            while t != n-1 && y != 1 && y != n-1
                y = (y*y)%n
                t *= 2
            end
            return false if y != n-1 && t & 1 == 0
        end
        return true
    end
end

module ModMath    # "http://ja.wikipedia.org/wiki/ミラー-ラビン素数判定法" より
    def ModMath.pow(base, power, mod)
        result = 1
        while power > 0
            result = (result*base)%mod if power & 1 == 1
            base = (base*base)%mod
            power >>= 1;
        end
        result
    end
end