#!/usr/local/bin/ruby -w

require 'benchmark' 

class Integer    
   def prime?         # cf. http://snippets.dzone.com/posts/show/4636
     n = self.abs()
     return true if n == 2
     return false if n == 1 || n & 1 == 0
     d = n-1
     d >>= 1 while d & 1 == 0
     20.times do                               # 20 = k from above
       a = rand(n-2) + 1
       t = d
       y = ModMath.pow(a,t,n)                  # implemented below
       while t != n-1 && y != 1 && y != n-1
         y = (y * y) % n
         t <<= 1
       end
       return false if y != n-1 && t & 1 == 0
     end
     return true
   end
end
 
module ModMath
   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



class Integer

   def primes   # cf. http://snippets.dzone.com/posts/show/3734

      sieve = []
      3.step(self, 2) { |i| sieve[i] = i }
      sieve[1] = nil
      sieve[2] = 2

      3.step(Math.sqrt(self).floor, 2) do |i| 
         next unless sieve[i]
         (i*i).step(self, i) do |j|
            sieve[j] = nil
         end
      end

      sieve.compact!

   end 


   def primes2       # cf. http://snippets.dzone.com/posts/show/3734

      primes = [nil, nil].concat((2..self).to_a)

      (2 .. Math.sqrt(self)).each do |i|
         next unless primes[i]
            (i*i).step(self, i) do |j|
               primes[j] = nil
            end
      end
	
      primes.compact!

   end

end



class Integer

   @primes_cache ||= 100.primes
   #@primes_cache ||= [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]
   class << self; attr_accessor :primes_cache; end

   @nthprime_limit ||= 5_000_000
   class << self; attr_reader :nthprime_limit; end
   #class << self; attr_accessor :nthprime_limit; end


   def sieve_size
      num = self.to_i.abs
      logn = Math.log(num)
      if num < 15985     # cf. http://primes.utm.edu/howmany.shtml
         ( num * (logn + Math.log(logn) - 1 + 1.8 * Math.log(logn) / logn) ).floor
      else
         ( num * (logn + Math.log(logn) - 0.9427) ).floor
      end
   end


   def nthprime

      num = self.to_i.abs

      # set a limit; cf. http://primes.utm.edu/nthprime/: The 5,000,000th prime is 86,028,121.
      raise "#{num}: number too big for Integer#nthprime" if num > Integer.nthprime_limit      

      primes_cache_size = Integer.primes_cache.size

      if num > primes_cache_size

         # optional; cf. Integer#nthprime_add and Integer#nthprime_add_mr below
         if num >= 50_000 && num < 100_000 && (num - primes_cache_size) < 5000 then return num.nthprime_add end
         if num >= 100_000 && num < 200_000 && (num - primes_cache_size) < 15000 then return num.nthprime_add end
         if num >= 200_000 && (num - primes_cache_size) < 35000 then return num.nthprime_add end
         
         logn = Math.log(num)

         if num < 15985       # cf. http://primes.utm.edu/howmany.shtml
            limit = ( num * (logn + Math.log(logn) - 1 + 1.8 * Math.log(logn) / logn) ).floor
            Integer.primes_cache = limit.primes
         else
            limit = ( num * (logn + Math.log(logn) - 0.9427) ).floor
            Integer.primes_cache = limit.primes
         end

=begin
         elsif num >= 15985 && num <= 1_000_000
            limit = ( num * (logn + Math.log(logn) - 0.9427) ).floor
            Integer.primes_cache = limit.primes
         elsif num > 1_000_000
            Integer.primes_cache = 20_000_000.primes
         end
=end

      else
         return Integer.primes_cache.at(num-1)
      end


      if num <= Integer.primes_cache.size
         return Integer.primes_cache.at(num-1)
      end

      if Integer.primes_cache.size > 2_500_000
         num.nthprime_add_mr 
      else
         num.nthprime_add   # faster for a (relatively) small prime cache
      end

   end


   def nthprime_add       #  add primes to Integer.primes_cache up to the nth prime

      num = self.to_i.abs

      if num <= Integer.primes_cache.size
         return Integer.primes_cache.at(num-1)
      end

      last_prime = Integer.primes_cache.last
      last_prime_divmod = last_prime.divmod(6)

      if last_prime_divmod.last == 1
         i = last_prime_divmod.first
         Integer.primes_cache.pop      # avoid a duplicate prime
      else
         i = last_prime_divmod.first + 1
      end


      while Integer.primes_cache.size < num

         n1 = 6*i+1       # cf. http://betterexplained.com/articles/another-look-at-prime-numbers/ and
         n2 = n1+4        # http://everything2.com/index.pl?node_id=1176369
         i += 1

         [n1, n2].each do |p| 

            next if p % 5 == 0 || p % 7 == 0

            next_num_sqrt = Math.sqrt(p).floor

            Integer.primes_cache.each do |d| 
               if d > next_num_sqrt 
                  Integer.primes_cache << p
                  #print "\r\e[0K#{Integer.primes_cache.size}"
                  #$stdout.flush
                  break
               elsif p % d == 0
                  break
               end  
            end
         end   # each  

      end  # while

      return Integer.primes_cache.at(num-1)

   end


   def nthprime_add_mr       #  add Miller-Rabin primes to Integer.primes_cache up to the nth prime

      num = self.to_i.abs

      if num <= Integer.primes_cache.size
         return Integer.primes_cache.at(num-1)
      end

      last_prime = Integer.primes_cache.last
      last_prime_divmod = last_prime.divmod(6)

      if last_prime_divmod.last == 1
         i = last_prime_divmod.first
         Integer.primes_cache.pop
      else
         i = last_prime_divmod.first + 1
      end


      while Integer.primes_cache.size < num

         search_next_prime = true

         while search_next_prime

            n1 = 6*i+1       
            n2 = n1+4        
            i += 1

            [n1, n2].each do |p| 

               next if p % 5 == 0 || p % 7 == 0

               if p.prime?
                  Integer.primes_cache << p
                  search_next_prime = false
                  #print "\r\e[0K#{Integer.primes_cache.size}"
                  #$stdout.flush
               end
            end
         end

      end  #  first while loop

      return Integer.primes_cache.at(num-1)

   end


#------------------------------------------- some additional prime methods


   def next_primes_in_cache    # next primes in Integer.primes_cache

      search_next_primes = self.to_i.abs

      last_prime = Integer.primes_cache.last
      last_prime_divmod = last_prime.divmod(6)

      if last_prime_divmod.last == 1
         i = last_prime_divmod.first
         Integer.primes_cache.pop
      else
         i = last_prime_divmod.first + 1
      end

      while search_next_primes > 0

         n1 = 6*i+1       
         n2 = n1+4        
         i += 1

         [n1, n2].each do |p| 

            next if p % 5 == 0 || p % 7 == 0
            next_num_sqrt = Math.sqrt(p).floor

            Integer.primes_cache.each do |d| 
               if d > next_num_sqrt 
                  Integer.primes_cache << p
                  search_next_primes -= 1
                  #print "\r\e[0K#{Integer.primes_cache.size}"
                  #$stdout.flush
                  break
               elsif p % d == 0
                  break
               end  
            end
         end 
      end   # while

      Integer.primes_cache.last(self.to_i.abs)

   end


   def next_mr_prime     # next Miller-Rabin prime
     
      last_num = self.to_i.abs
      last_num_divmod = last_num.divmod(6)

      if last_num_divmod.last == 1
         i = last_num_divmod.first
      else
         i = last_num_divmod.first + 1
      end

      next_prime = nil
      search_next_prime = true

      while search_next_prime

         n1 = 6*i+1       
         n2 = n1+4        
         i += 1

         [n1, n2].each do |p| 

            next if p % 5 == 0 || p % 7 == 0

            if p > last_num && p.prime?
               next_prime = p
               search_next_prime = false 
               break
            end
         end
      end

         next_prime

   end



   def next_mr_primes(n)     # next Miller-Rabin primes

      last_num = self.to_i.abs
      last_num_divmod = last_num.divmod(6)

      if last_num_divmod.last == 1
         i = last_num_divmod.first
      else
         i = last_num_divmod.first + 1
      end

      next_primes = []
      search_next_primes = n.to_i.abs

      while search_next_primes > 0

         n1 = 6*i+1      
         n2 = n1+4       
         i += 1

         [n1, n2].each do |p| 

            next if p % 5 == 0 || p % 7 == 0

            if p > last_num && p.prime?
               next_primes << p
               search_next_primes -= 1 
            end
         end
      end

      next_primes

   end

end



num1 = 10_000
num1 = 1_000
num1 = 5_000

num2 = 210_349
num2 = 100_125
num2 = 55_000

ret1 = nil
ret2 = nil
ret3 = nil
ret4 = nil


Benchmark.bm(16) do |t| 

   t.report("first #{num1} primes: ") do
      ret1 = num1.nthprime_add
   end 

   Integer.primes_cache.clear
   Integer.primes_cache.concat(100.primes)

   t.report("first #{num1} primes: ") do
      ret2 = num1.nthprime_add_mr
   end 

   Integer.primes_cache.clear
   Integer.primes_cache.concat(100.primes)

   t.report("first #{num1} primes: ") do
      ret3 = num1.nthprime
   end 

   Integer.primes_cache.clear
   Integer.primes_cache.concat(100.primes)

   t.report("first #{num2} primes: ") do
      ret4 = num2.nthprime
   end 

end


puts
puts "the #{num1}th prime number: #{ret1}"
puts "the #{num1}th prime number: #{ret2}"
puts "the #{num1}th prime number: #{ret3}"
puts "the #{num2}th prime number: #{ret4}"
puts
puts "the #{num1}th prime: #{Integer.primes_cache.at(num1-1)}"
puts "the #{num2}th prime: #{Integer.primes_cache.at(num2-1)}"
puts
p Integer.primes_cache.first(10)
p Integer.primes_cache.last(10)
p Integer.primes_cache.size
puts


p 15.next_primes_in_cache

puts 594_213.next_mr_prime

p 149_137.next_mr_primes(5)

puts
p 1_000.sieve_size
p 1_000_000.sieve_size
p 2_500_000.sieve_size
p 5_000_000.sieve_size
p 10_000_000.sieve_size
p 100_000_000.sieve_size


=begin

# check with primegen.c, http://cr.yp.to/primegen.html
primes 2 48611 | nl | tail -n 1
primes 2 104729 | nl | tail -n 1
primes 2 679277 | nl | tail -n 1
primes 2 1301423 | nl | tail -n 1
primes 2 2903603 | nl | tail -n 1
primes 1303129 1303283 | nl
primes 594213 $((594213 + 500)) | head -n 1
primes 149137 $((149137 + 1000)) | head -n 5

# further prime tools from http://libtom.org
- LibTomMath
   bn_mp_prime_is_prime.c
   bn_mp_prime_miller_rabin.c
- TomsFastMath
   fp_isprime.c
   fp_prime_miller_rabin.c

=end

# From: http://en.wikipedia.org/wiki/Miller-Rabin_primality_test

class Integer
   def prime?
     n = self.abs()
     return true if n == 2
     return false if n == 1 || n & 1 == 0

     # cf. http://betterexplained.com/articles/another-look-at-prime-numbers/ and
     # http://everything2.com/index.pl?node_id=1176369

     return false if n > 3 && n % 6 != 1 && n % 6 != 5     # added

     d = n-1
     d >>= 1 while d & 1 == 0
     20.times do                               # 20 = k from above
       a = rand(n-2) + 1
       t = d
       y = ModMath.pow(a,t,n)                  # implemented below
       while t != n-1 && y != 1 && y != n-1
         y = (y * y) % n
         t <<= 1
       end
       return false if y != n-1 && t & 1 == 0
     end
     return true
   end
end
 
module ModMath
   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


#------------------------------------------------------------------


# From: http://www.h2np.net/tips/ruby-cipher-math.txt
# Author: Hironobu SUZUKI

# $Id: mathh.rb,v 1.5 2002/11/24 16:15:56 hironobu Exp $
# Simple Secure Stream Cryptography
# Mathematical Library
# Hironobu SUZUKI <hironobu@h2np.net>
# LGPL, (C) 2002, Hironobu SUZUKI

class Random  
  def initialize()
    @random_device="/dev/urandom"
    @verbose=false
  end
  def random_device=(str)
    @random_device=str
  end
  def verbose
    @verbose=true
  end
  def get(size)
    if (@verbose) ;STDERR.print "+";end
    value=0
    r=File.open(@random_device)
    r.read(size).each_byte {|c|
      value = (value << 8) + c.to_i
    }
    r.close
    return value
    
  end
  def get_prime(lower,higher)
    bytesize = higher.bytesize
    while true
      r = self.get(bytesize)
      if (r % 2 == 0) ; r -= 1 ;   end
      if ( lower <=  r && r <= higher && r.prime? == true )
	return r
      end
    end
  end
end


class Integer
  def cdiv(d)			# ceil divide
    w=self.divmod(d)
    if w[1] > 0 ;  w[0] +=1 ; end
    return w[0]
  end
  def powm(i,n)	# square-and-multiply for exp modulo n
    if ( i == 0 )
      return 1
    end
    b=1
    a=self

    if ( (i & 1) == 1 )
      b=a
    end
    i = i>>1

    while ( i > 0 )
      a = (a * a) % n
      if ( (i & 1) == 1 )
	b = (a * b) % n 
      end
      i = i>>1
    end
    return b
  end
  def gcd(b)			# GCD
    a = self
    if ( a < b ); t=b; b=a; a=t; end
    while b > 0 
      r = a % b 
      a = b
      b = r
    end
    return a
  end
  def invert(n)  # Extension Euclid Algorithm
    a = n        # See Knuth's The Art of Computer Programming 
    b = self % n # Vol2 pp.342 -- 343  
    p = 0 ; q = 1;  v = 0;  u = 1;
    while q > 0 
      p = a / b
      q =  a % b
      w = v - (p*u)
      ## DEBUG    printf "%d = %d(%d) + %d  (%d)\n", a,p,b,q,v ###
      a = b
      b = q
      v = u
      u = w
    end
    return (v+n) % n
  end

# Miller-Rabin Test  (Prime Test)
# See, http://www.cs.albany.edu/~berg/csi445/Assignments/pa4.html
  def bitarray(n) 
    b=Array::new  
    i=0       
    v=n
    while v > 0
      b[i] = (0x1 & v)
      v = v/2
      i = i + 1
    end
    return b
  end
  def miller_rabin(n,s)
    b=bitarray(n-1)
    i=b.size 
    j =1
    while j <= s
      a = 1 + (rand(n-2).to_i)
      if witness(a,n,i,b) == true
	return false
      end
      j+=1
    end
    return true
  end
  def witness(a,n,i,b)
    d=1
    x=0
    while i > 0 
      x = d 
      d = (d**2) % n
      if ( (d == 1) && (x != 1) && (x != (n-1)) )
	return true
      end
      if ( b[i-1] == 1 )
	d = (d * a ) % n
      end
      i -= 1
    end
    if ( d != 1) 
      return true
    end
    return false
  end

  #def prime?
  def is_prime?
    a = self
    return miller_rabin(a,30)    
  end

  def bytesize
    n = self
    i=0
    while n > 0
      n = n >> 8
      i += 1
    end
    return i
  end
  def bitsize
    n = self
    i=0
    while n > 0
      n = n >> 1
      i += 1
    end
    return i
  end
end


p 107565456790871.prime?      # true
p 107565456790871.is_prime?   # true

a = 107565456790871

begin
  a += 2
end while !a.prime? && !a.is_prime?

puts a   #=> 107565456790991 (next prime)

CharExponent = 3
BitsPerChar = 1 << CharExponent
LowMask = BitsPerChar - 1

def sieve(m)
  ret = []
  items = "\xFF" * ((m / BitsPerChar) + 1)
  masks = ""
  BitsPerChar.times do |b|
    masks << (1 << b).chr
  end

  count = 0
  pmax = m - 1
  2.step(pmax, 1) do |p|
    if items[p >> CharExponent][p & LowMask] == 1
      ret << p
      count += 1
      p.step(pmax, p) do |mult|
	a = mult >> CharExponent
	b = mult & LowMask
	items[a] -= masks[b] if items[a][b] != 0
        #p items
      end
    end
  end
  #count
  ret
end

n = (ARGV[0] || 2).to_i
n.step(n - 2, -1) do |exponent|
  break if exponent < 0
  m = 2 ** exponent * 10_000
  primes = sieve(m)
  count = primes.size
  puts
  printf "Primes up to %8d %8d\n", m, count
  puts "last ten primes: #{primes[-10..-1].join(' ')}"
  puts
end

import sys
import random

def toBinary(n):
  r = []
  while (n > 0):
    r.append(n % 2)
    n = n / 2
  return r

def test(a, n):
  """
  test(a, n) -> bool Tests whether n is complex.

  Returns:
    - True, if n is complex.
    - False, if n is probably prime.
  """
  b = toBinary(n - 1)
  d = 1
  for i in xrange(len(b) - 1, -1, -1):
    x = d
    d = (d * d) % n
    if d == 1 and x != 1 and x != n - 1:
      return True # Complex
    if b[i] == 1:
      d = (d * a) % n
  if d != 1:
    return True # Complex
  return False # Prime

def MillerRabin(n, s = 50):
  """
    MillerRabin(n, s = 1000) -> bool Checks whether n is prime or not

    Returns:
      - True, if n is probably prime.
      - False, if n is complex.
  """
  for j in xrange(1, s + 1):
    a = random.randint(1, n - 1)
    if (test(a, n)):
      return False # n is complex
  return True # n is prime

def main(argv):
  print MillerRabin(int(argv[0]), int(argv[1]))

if __name__ == "__main__":
  main(sys.argv[1:])

def euclidExtended(a, b):
  if b == 0:
    return a, 1, 0
  dd, xx, yy = euclidExtended(b, a % b)
  d, x, y = dd, yy, xx - int(a / b) * yy
  return d, x, y

int phi(int n) {
  double rez = n;

  int i = 0;
  while ((i < nrPrime) && (prime[i] <= n)) {
    if (n % prime[i] == 0) {
      rez *= (double)(prime[i] - 1) / (double)prime[i];
    }
    ++i;
  }

  return rez;
}

#include <string.h>

int prime[78498]; // 78498 is the number of prime numbers smaller than 1 000 000.
int nrPrime = 0;

void genPrime(int n) {
  char p[1000001];

  memset(p, 1, sizeof(char) * (n + 1));

  p[1] = 0;

  int i = 2;
  for (; i <= n; ++i)
    if (p[i]) {
      prime[nrPrime] = i;
      ++nrPrime;
      int j = i * 2;
      while (j <= n) {
        p[j] = 0;
        j += i;
      }
    }
}

max = Integer(ARGV.shift || 100)
	
sieve = [nil, nil] + (2 .. max).to_a
	
(2 .. Math.sqrt(max)).each do |i|
  next unless sieve[i]
  (i*i).step(max, i) do |j|
    sieve[j] = nil
  end
end
	
puts sieve.compact.join(", ")



#------------------------------------



#!/usr/local/bin/ruby -w

require 'benchmark' 

class Integer

   def primes1

      primes = [nil, nil].concat((2..self).to_a)

      (2 .. Math.sqrt(self).floor).each do |i|
         next unless primes[i]
            (i*i).step(self, i) do |j|
               primes[j] = nil
            end
      end
	
      primes.compact!

   end  


   def primes2

      sieve = []
      3.step(self, 2) { |i| sieve[i] = i }
      sieve[1] = nil
      sieve[2] = 2

      3.step(Math.sqrt(self).floor, 2) do |i| 
         next unless sieve[i]
         (i*i).step(self, i) do |j|
            sieve[j] = nil
         end
      end

      sieve.compact!

   end 


   def primes3

      n2 = 0
      i = 0
      sieve = []

      while n2 < self
         n1 = 6*i+1       # cf. http://betterexplained.com/articles/another-look-at-prime-numbers/ and
         n2 = n1+4        # http://everything2.com/index.pl?node_id=1176369
         i += 1

         sieve[n1] = n1
         sieve[n2] = n2

      end  

      sieve[1] = nil
      sieve[2] = 2
      sieve[3] = 3
      sieve[-1] = nil

      5.step(Math.sqrt(self).floor, 2) do |i| 
         next unless sieve[i]
         (i*i).step(self, i) do |j|     
            sieve[j] = nil
         end