1  
   2  
   3  #!/usr/local/bin/ruby -w
   4  
   5  require 'benchmark' 
   6  
   7  class Integer    
   8     def prime?         # cf. http://snippets.dzone.com/posts/show/4636
   9       n = self.abs()
  10       return true if n == 2
  11       return false if n == 1 || n & 1 == 0
  12       d = n-1
  13       d >>= 1 while d & 1 == 0
  14       20.times do                               # 20 = k from above
  15         a = rand(n-2) + 1
  16         t = d
  17         y = ModMath.pow(a,t,n)                  # implemented below
  18         while t != n-1 && y != 1 && y != n-1
  19           y = (y * y) % n
  20           t <<= 1
  21         end
  22         return false if y != n-1 && t & 1 == 0
  23       end
  24       return true
  25     end
  26  end
  27   
  28  module ModMath
  29     def ModMath.pow(base, power, mod)
  30       result = 1
  31       while power > 0
  32         result = (result * base) % mod if power & 1 == 1
  33         base = (base * base) % mod
  34         power >>= 1;
  35       end
  36       result
  37     end
  38  end
  39  
  40  
  41  
  42  class Integer
  43  
  44     def primes   # cf. http://snippets.dzone.com/posts/show/3734
  45  
  46        sieve = []
  47        3.step(self, 2) { |i| sieve[i] = i }
  48        sieve[1] = nil
  49        sieve[2] = 2
  50  
  51        3.step(Math.sqrt(self).floor, 2) do |i| 
  52           next unless sieve[i]
  53           (i*i).step(self, i) do |j|
  54              sieve[j] = nil
  55           end
  56        end
  57  
  58        sieve.compact!
  59  
  60     end 
  61  
  62  
  63     def primes2       # cf. http://snippets.dzone.com/posts/show/3734
  64  
  65        primes = [nil, nil].concat((2..self).to_a)
  66  
  67        (2 .. Math.sqrt(self)).each do |i|
  68           next unless primes[i]
  69              (i*i).step(self, i) do |j|
  70                 primes[j] = nil
  71              end
  72        end
  73  	
  74        primes.compact!
  75  
  76     end
  77  
  78  end
  79  
  80  
  81  
  82  class Integer
  83  
  84     @primes_cache ||= 100.primes
  85     #@primes_cache ||= [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]
  86     class << self; attr_accessor :primes_cache; end
  87  
  88     @nthprime_limit ||= 5_000_000
  89     class << self; attr_reader :nthprime_limit; end
  90     #class << self; attr_accessor :nthprime_limit; end
  91  
  92  
  93     def sieve_size
  94        num = self.to_i.abs
  95        logn = Math.log(num)
  96        if num < 15985     # cf. http://primes.utm.edu/howmany.shtml
  97           ( num * (logn + Math.log(logn) - 1 + 1.8 * Math.log(logn) / logn) ).floor
  98        else
  99           ( num * (logn + Math.log(logn) - 0.9427) ).floor
 100        end
 101     end
 102  
 103  
 104     def nthprime
 105  
 106        num = self.to_i.abs
 107  
 108        # set a limit; cf. http://primes.utm.edu/nthprime/: The 5,000,000th prime is 86,028,121.
 109        raise "#{num}: number too big for Integer#nthprime" if num > Integer.nthprime_limit      
 110  
 111        primes_cache_size = Integer.primes_cache.size
 112  
 113        if num > primes_cache_size
 114  
 115           # optional; cf. Integer#nthprime_add and Integer#nthprime_add_mr below
 116           if num >= 50_000 && num < 100_000 && (num - primes_cache_size) < 5000 then return num.nthprime_add end
 117           if num >= 100_000 && num < 200_000 && (num - primes_cache_size) < 15000 then return num.nthprime_add end
 118           if num >= 200_000 && (num - primes_cache_size) < 35000 then return num.nthprime_add end
 119           
 120           logn = Math.log(num)
 121  
 122           if num < 15985       # cf. http://primes.utm.edu/howmany.shtml
 123              limit = ( num * (logn + Math.log(logn) - 1 + 1.8 * Math.log(logn) / logn) ).floor
 124              Integer.primes_cache = limit.primes
 125           else
 126              limit = ( num * (logn + Math.log(logn) - 0.9427) ).floor
 127              Integer.primes_cache = limit.primes
 128           end
 129  
 130  =begin
 131           elsif num >= 15985 && num <= 1_000_000
 132              limit = ( num * (logn + Math.log(logn) - 0.9427) ).floor
 133              Integer.primes_cache = limit.primes
 134           elsif num > 1_000_000
 135              Integer.primes_cache = 20_000_000.primes
 136           end
 137  =end
 138  
 139        else
 140           return Integer.primes_cache.at(num-1)
 141        end
 142  
 143  
 144        if num <= Integer.primes_cache.size
 145           return Integer.primes_cache.at(num-1)
 146        end
 147  
 148        if Integer.primes_cache.size > 2_500_000
 149           num.nthprime_add_mr 
 150        else
 151           num.nthprime_add   # faster for a (relatively) small prime cache
 152        end
 153  
 154     end
 155  
 156  
 157     def nthprime_add       #  add primes to Integer.primes_cache up to the nth prime
 158  
 159        num = self.to_i.abs
 160  
 161        if num <= Integer.primes_cache.size
 162           return Integer.primes_cache.at(num-1)
 163        end
 164  
 165        last_prime = Integer.primes_cache.last
 166        last_prime_divmod = last_prime.divmod(6)
 167  
 168        if last_prime_divmod.last == 1
 169           i = last_prime_divmod.first
 170           Integer.primes_cache.pop      # avoid a duplicate prime
 171        else
 172           i = last_prime_divmod.first + 1
 173        end
 174  
 175  
 176        while Integer.primes_cache.size < num
 177  
 178           n1 = 6*i+1       # cf. http://betterexplained.com/articles/another-look-at-prime-numbers/ and
 179           n2 = n1+4        # http://everything2.com/index.pl?node_id=1176369
 180           i += 1
 181  
 182           [n1, n2].each do |p| 
 183  
 184              next if p % 5 == 0 || p % 7 == 0
 185  
 186              next_num_sqrt = Math.sqrt(p).floor
 187  
 188              Integer.primes_cache.each do |d| 
 189                 if d > next_num_sqrt 
 190                    Integer.primes_cache << p
 191                    #print "\r\e[0K#{Integer.primes_cache.size}"
 192                    #$stdout.flush
 193                    break
 194                 elsif p % d == 0
 195                    break
 196                 end  
 197              end
 198           end   # each  
 199  
 200        end  # while
 201  
 202        return Integer.primes_cache.at(num-1)
 203  
 204     end
 205  
 206  
 207     def nthprime_add_mr       #  add Miller-Rabin primes to Integer.primes_cache up to the nth prime
 208  
 209        num = self.to_i.abs
 210  
 211        if num <= Integer.primes_cache.size
 212           return Integer.primes_cache.at(num-1)
 213        end
 214  
 215        last_prime = Integer.primes_cache.last
 216        last_prime_divmod = last_prime.divmod(6)
 217  
 218        if last_prime_divmod.last == 1
 219           i = last_prime_divmod.first
 220           Integer.primes_cache.pop
 221        else
 222           i = last_prime_divmod.first + 1
 223        end
 224  
 225  
 226        while Integer.primes_cache.size < num
 227  
 228           search_next_prime = true
 229  
 230           while search_next_prime
 231  
 232              n1 = 6*i+1       
 233              n2 = n1+4        
 234              i += 1
 235  
 236              [n1, n2].each do |p| 
 237  
 238                 next if p % 5 == 0 || p % 7 == 0
 239  
 240                 if p.prime?
 241                    Integer.primes_cache << p
 242                    search_next_prime = false
 243                    #print "\r\e[0K#{Integer.primes_cache.size}"
 244                    #$stdout.flush
 245                 end
 246              end
 247           end
 248  
 249        end  #  first while loop
 250  
 251        return Integer.primes_cache.at(num-1)
 252  
 253     end
 254  
 255  
 256  #------------------------------------------- some additional prime methods
 257  
 258  
 259     def next_primes_in_cache    # next primes in Integer.primes_cache
 260  
 261        search_next_primes = self.to_i.abs
 262  
 263        last_prime = Integer.primes_cache.last
 264        last_prime_divmod = last_prime.divmod(6)
 265  
 266        if last_prime_divmod.last == 1
 267           i = last_prime_divmod.first
 268           Integer.primes_cache.pop
 269        else
 270           i = last_prime_divmod.first + 1
 271        end
 272  
 273        while search_next_primes > 0
 274  
 275           n1 = 6*i+1       
 276           n2 = n1+4        
 277           i += 1
 278  
 279           [n1, n2].each do |p| 
 280  
 281              next if p % 5 == 0 || p % 7 == 0
 282              next_num_sqrt = Math.sqrt(p).floor
 283  
 284              Integer.primes_cache.each do |d| 
 285                 if d > next_num_sqrt 
 286                    Integer.primes_cache << p
 287                    search_next_primes -= 1
 288                    #print "\r\e[0K#{Integer.primes_cache.size}"
 289                    #$stdout.flush
 290                    break
 291                 elsif p % d == 0
 292                    break
 293                 end  
 294              end
 295           end 
 296        end   # while
 297  
 298        Integer.primes_cache.last(self.to_i.abs)
 299  
 300     end
 301  
 302  
 303     def next_mr_prime     # next Miller-Rabin prime
 304       
 305        last_num = self.to_i.abs
 306        last_num_divmod = last_num.divmod(6)
 307  
 308        if last_num_divmod.last == 1
 309           i = last_num_divmod.first
 310        else
 311           i = last_num_divmod.first + 1
 312        end
 313  
 314        next_prime = nil
 315        search_next_prime = true
 316  
 317        while search_next_prime
 318  
 319           n1 = 6*i+1       
 320           n2 = n1+4        
 321           i += 1
 322  
 323           [n1, n2].each do |p| 
 324  
 325              next if p % 5 == 0 || p % 7 == 0
 326  
 327              if p > last_num && p.prime?
 328                 next_prime = p
 329                 search_next_prime = false 
 330                 break
 331              end
 332           end
 333        end
 334  
 335           next_prime
 336  
 337     end
 338  
 339  
 340  
 341     def next_mr_primes(n)     # next Miller-Rabin primes
 342  
 343        last_num = self.to_i.abs
 344        last_num_divmod = last_num.divmod(6)
 345  
 346        if last_num_divmod.last == 1
 347           i = last_num_divmod.first
 348        else
 349           i = last_num_divmod.first + 1
 350        end
 351  
 352        next_primes = []
 353        search_next_primes = n.to_i.abs
 354  
 355        while search_next_primes > 0
 356  
 357           n1 = 6*i+1      
 358           n2 = n1+4       
 359           i += 1
 360  
 361           [n1, n2].each do |p| 
 362  
 363              next if p % 5 == 0 || p % 7 == 0
 364  
 365              if p > last_num && p.prime?
 366                 next_primes << p
 367                 search_next_primes -= 1 
 368              end
 369           end
 370        end
 371  
 372        next_primes
 373  
 374     end
 375  
 376  end
 377  
 378  
 379  
 380  num1 = 10_000
 381  num1 = 1_000
 382  num1 = 5_000
 383  
 384  num2 = 210_349
 385  num2 = 100_125
 386  num2 = 55_000
 387  
 388  ret1 = nil
 389  ret2 = nil
 390  ret3 = nil
 391  ret4 = nil
 392  
 393  
 394  Benchmark.bm(16) do |t| 
 395  
 396     t.report("first #{num1} primes: ") do
 397        ret1 = num1.nthprime_add
 398     end 
 399  
 400     Integer.primes_cache.clear
 401     Integer.primes_cache.concat(100.primes)
 402  
 403     t.report("first #{num1} primes: ") do
 404        ret2 = num1.nthprime_add_mr
 405     end 
 406  
 407     Integer.primes_cache.clear
 408     Integer.primes_cache.concat(100.primes)
 409  
 410     t.report("first #{num1} primes: ") do
 411        ret3 = num1.nthprime
 412     end 
 413  
 414     Integer.primes_cache.clear
 415     Integer.primes_cache.concat(100.primes)
 416  
 417     t.report("first #{num2} primes: ") do
 418        ret4 = num2.nthprime
 419     end 
 420  
 421  end
 422  
 423  
 424  puts
 425  puts "the #{num1}th prime number: #{ret1}"
 426  puts "the #{num1}th prime number: #{ret2}"
 427  puts "the #{num1}th prime number: #{ret3}"
 428  puts "the #{num2}th prime number: #{ret4}"
 429  puts
 430  puts "the #{num1}th prime: #{Integer.primes_cache.at(num1-1)}"
 431  puts "the #{num2}th prime: #{Integer.primes_cache.at(num2-1)}"
 432  puts
 433  p Integer.primes_cache.first(10)
 434  p Integer.primes_cache.last(10)
 435  p Integer.primes_cache.size
 436  puts
 437  
 438  
 439  p 15.next_primes_in_cache
 440  
 441  puts 594_213.next_mr_prime
 442  
 443  p 149_137.next_mr_primes(5)
 444  
 445  puts
 446  p 1_000.sieve_size
 447  p 1_000_000.sieve_size
 448  p 2_500_000.sieve_size
 449  p 5_000_000.sieve_size
 450  p 10_000_000.sieve_size
 451  p 100_000_000.sieve_size
 452  
 453  
 454  =begin
 455  
 456  # check with primegen.c, http://cr.yp.to/primegen.html
 457  primes 2 48611 | nl | tail -n 1
 458  primes 2 104729 | nl | tail -n 1
 459  primes 2 679277 | nl | tail -n 1
 460  primes 2 1301423 | nl | tail -n 1
 461  primes 2 2903603 | nl | tail -n 1
 462  primes 1303129 1303283 | nl
 463  primes 594213 $((594213 + 500)) | head -n 1
 464  primes 149137 $((149137 + 1000)) | head -n 5
 465  
 466  # further prime tools from http://libtom.org
 467  - LibTomMath
 468     bn_mp_prime_is_prime.c
 469     bn_mp_prime_miller_rabin.c
 470  - TomsFastMath
 471     fp_isprime.c
 472     fp_prime_miller_rabin.c
 473  
 474  =end

   1  
   2  
   3  # From: http://en.wikipedia.org/wiki/Miller-Rabin_primality_test
   4  
   5  class Integer
   6     def prime?
   7       n = self.abs()
   8       return true if n == 2
   9       return false if n == 1 || n & 1 == 0
  10  
  11       # cf. http://betterexplained.com/articles/another-look-at-prime-numbers/ and
  12       # http://everything2.com/index.pl?node_id=1176369
  13  
  14       return false if n > 3 && n % 6 != 1 && n % 6 != 5     # added
  15  
  16       d = n-1
  17       d >>= 1 while d & 1 == 0
  18       20.times do                               # 20 = k from above
  19         a = rand(n-2) + 1
  20         t = d
  21         y = ModMath.pow(a,t,n)                  # implemented below
  22         while t != n-1 && y != 1 && y != n-1
  23           y = (y * y) % n
  24           t <<= 1
  25         end
  26         return false if y != n-1 && t & 1 == 0
  27       end
  28       return true
  29     end
  30  end
  31   
  32  module ModMath
  33     def ModMath.pow(base, power, mod)
  34       result = 1
  35       while power > 0
  36         result = (result * base) % mod if power & 1 == 1
  37         base = (base * base) % mod
  38         power >>= 1;
  39       end
  40       result
  41     end
  42  end
  43  
  44  
  45  #------------------------------------------------------------------
  46  
  47  
  48  # From: http://www.h2np.net/tips/ruby-cipher-math.txt
  49  # Author: Hironobu SUZUKI
  50  
  51  # $Id: mathh.rb,v 1.5 2002/11/24 16:15:56 hironobu Exp $
  52  # Simple Secure Stream Cryptography
  53  # Mathematical Library
  54  # Hironobu SUZUKI <hironobu@h2np.net>
  55  # LGPL, (C) 2002, Hironobu SUZUKI
  56  
  57  class Random  
  58    def initialize()
  59      @random_device="/dev/urandom"
  60      @verbose=false
  61    end
  62    def random_device=(str)
  63      @random_device=str
  64    end
  65    def verbose
  66      @verbose=true
  67    end
  68    def get(size)
  69      if (@verbose) ;STDERR.print "+";end
  70      value=0
  71      r=File.open(@random_device)
  72      r.read(size).each_byte {|c|
  73        value = (value << 8) + c.to_i
  74      }
  75      r.close
  76      return value
  77      
  78    end
  79    def get_prime(lower,higher)
  80      bytesize = higher.bytesize
  81      while true
  82        r = self.get(bytesize)
  83        if (r % 2 == 0) ; r -= 1 ;   end
  84        if ( lower <=  r && r <= higher && r.prime? == true )
  85  	return r
  86        end
  87      end
  88    end
  89  end
  90  
  91  
  92  class Integer
  93    def cdiv(d)			# ceil divide
  94      w=self.divmod(d)
  95      if w[1] > 0 ;  w[0] +=1 ; end
  96      return w[0]
  97    end
  98    def powm(i,n)	# square-and-multiply for exp modulo n
  99      if ( i == 0 )
 100        return 1
 101      end
 102      b=1
 103      a=self
 104  
 105      if ( (i & 1) == 1 )
 106        b=a
 107      end
 108      i = i>>1
 109  
 110      while ( i > 0 )
 111        a = (a * a) % n
 112        if ( (i & 1) == 1 )
 113  	b = (a * b) % n 
 114        end
 115        i = i>>1
 116      end
 117      return b
 118    end
 119    def gcd(b)			# GCD
 120      a = self
 121      if ( a < b ); t=b; b=a; a=t; end
 122      while b > 0 
 123        r = a % b 
 124        a = b
 125        b = r
 126      end
 127      return a
 128    end
 129    def invert(n)  # Extension Euclid Algorithm
 130      a = n        # See Knuth's The Art of Computer Programming 
 131      b = self % n # Vol2 pp.342 -- 343  
 132      p = 0 ; q = 1;  v = 0;  u = 1;
 133      while q > 0 
 134        p = a / b
 135        q =  a % b
 136        w = v - (p*u)
 137        ## DEBUG    printf "%d = %d(%d) + %d  (%d)\n", a,p,b,q,v ###
 138        a = b
 139        b = q
 140        v = u
 141        u = w
 142      end
 143      return (v+n) % n
 144    end
 145  
 146  # Miller-Rabin Test  (Prime Test)
 147  # See, http://www.cs.albany.edu/~berg/csi445/Assignments/pa4.html
 148    def bitarray(n) 
 149      b=Array::new  
 150      i=0       
 151      v=n
 152      while v > 0
 153        b[i] = (0x1 & v)
 154        v = v/2
 155        i = i + 1
 156      end
 157      return b
 158    end
 159    def miller_rabin(n,s)
 160      b=bitarray(n-1)
 161      i=b.size 
 162      j =1
 163      while j <= s
 164        a = 1 + (rand(n-2).to_i)
 165        if witness(a,n,i,b) == true
 166  	return false
 167        end
 168        j+=1
 169      end
 170      return true
 171    end
 172    def witness(a,n,i,b)
 173      d=1
 174      x=0
 175      while i > 0 
 176        x = d 
 177        d = (d**2) % n
 178        if ( (d == 1) && (x != 1) && (x != (n-1)) )
 179  	return true
 180        end
 181        if ( b[i-1] == 1 )
 182  	d = (d * a ) % n
 183        end
 184        i -= 1
 185      end
 186      if ( d != 1) 
 187        return true
 188      end
 189      return false
 190    end
 191  
 192    #def prime?
 193    def is_prime?
 194      a = self
 195      return miller_rabin(a,30)    
 196    end
 197  
 198    def bytesize
 199      n = self
 200      i=0
 201      while n > 0
 202        n = n >> 8
 203        i += 1
 204      end
 205      return i
 206    end
 207    def bitsize
 208      n = self
 209      i=0
 210      while n > 0
 211        n = n >> 1
 212        i += 1
 213      end
 214      return i
 215    end
 216  end
 217  
 218  
 219  p 107565456790871.prime?      # true
 220  p 107565456790871.is_prime?   # true
 221  
 222  a = 107565456790871
 223  
 224  begin
 225    a += 2
 226  end while !a.prime? && !a.is_prime?
 227  
 228  puts a   #=> 107565456790991 (next prime)
 229  

   1  
   2  
   3  CharExponent = 3
   4  BitsPerChar = 1 << CharExponent
   5  LowMask = BitsPerChar - 1
   6  
   7  def sieve(m)
   8    ret = []
   9    items = "\xFF" * ((m / BitsPerChar) + 1)
  10    masks = ""
  11    BitsPerChar.times do |b|
  12      masks << (1 << b).chr
  13    end
  14  
  15    count = 0
  16    pmax = m - 1
  17    2.step(pmax, 1) do |p|
  18      if items[p >> CharExponent][p & LowMask] == 1
  19        ret << p
  20        count += 1
  21        p.step(pmax, p) do |mult|
  22  	a = mult >> CharExponent
  23  	b = mult & LowMask
  24  	items[a] -= masks[b] if items[a][b] != 0
  25          #p items
  26        end
  27      end
  28    end
  29    #count
  30    ret
  31  end
  32  
  33  n = (ARGV[0] || 2).to_i
  34  n.step(n - 2, -1) do |exponent|
  35    break if exponent < 0
  36    m = 2 ** exponent * 10_000
  37    primes = sieve(m)
  38    count = primes.size
  39    puts
  40    printf "Primes up to %8d %8d\n", m, count
  41    puts "last ten primes: #{primes[-10..-1].join(' ')}"
  42    puts
  43  end
  44  

   1  
   2  
   3  max = Integer(ARGV.shift || 100)
   4  	
   5  sieve = [nil, nil] + (2 .. max).to_a
   6  	
   7  (2 .. Math.sqrt(max)).each do |i|
   8    next unless sieve[i]
   9    (i*i).step(max, i) do |j|
  10      sieve[j] = nil
  11    end
  12  end
  13  	
  14  puts sieve.compact.join(", ")
  15  
  16  
  17  
  18  #------------------------------------
  19  
  20  
  21  
  22  #!/usr/local/bin/ruby -w
  23  
  24  require 'benchmark' 
  25  
  26  class Integer
  27  
  28     def primes1
  29  
  30        primes = [nil, nil].